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ABSTRACT 


In  many  strategic  shallow  water  areas  the  geoacoustic  properties  of  the 
sub-bottom  are  largely  unknown.  In  this  thesis  it  is  demonstrated  that 
inverse  theory  and  measured  data  from  a  single  hydrophone  can  be  used  to 
accurately  deduce  the  geoacoustic  properties  of  the  sub-bottom,  even  when 
the  initial  background  geoacoustic  model  is  a  highly  inaccurate  “guess”. 
Since  propagation  in  shallow  water  is  very  sensitive  to  the  geoacoustic 
properties  of  the  sub-bottom,  the  inverse  technique  developed  in  this  thesis 
presents  the  Navy  with  a  vitally  important,  practical,  and  inexpensive  means 
to  improve  sonar  performance  prediction  in  a  potentially  hostile  environment. 

To  provide  ground  truth  for  the  inverse  technique,  measured  data 
collected  during  Project  GEMINI  were  compared  to  the  inverse  solutions. 
Detailed,  site-specific  geoacoustic  models  were  developed  for  two  array 
locations  and  the  Finite  Element  Parabolic  Equation  (FEPE)  model  was  used 
to  estimate  transmission  loss  (TL).  The  model  estimates  firom  FEPE 
compared  well  with  the  measured  data  and  the  detailed  geoacoustic  models 
were  considered  as  “groimd  truth”.  To  test  the  efficacy  of  the  technique, 
initial  background  geoacoustic  models  were  constructed  assuming  no  a  priori 
information  of  the  bottom.  The  resultant  inverse  solution  was  used  to  predict 
the  geoacoustic  properties  at  each  of  the  sites.  The  final  results  were  in 
excellent  agreement  with  the  measured  data  and  the  resulting  inverse 
technique  TL  estimates  were  as  good  or  better  than  the  TL  estimates 
obtained  from  the  detailed,  site-specific  geoacoustic  models. 
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1.  INTRODUCTION 


A.  NEW  PRIORITIES 

The  United  States  Navy  (USN)  strategic  vision  articulated  in  "...From 
the  Sea, "(O’Keefe  et  al.,  1992)  asserts  U.S.  naval  forces  as  full  participants  in 
the  new  national  strategy  for  the  post-Cold  War  Era.  No  longer  focused 
primarily  on  open  ocean  global  conflict,  the  navsil  service’s  new  emphasis  is 
clearly  placed  on  regional  contingencies  in  the  littoral  regions  of  the  world. 
Accordingly,  the  Navy  has  downgraded  its  previous  anti-submarine  warfare 
(ASW)  requirements  that  stressed  blue-water  capabihties  against  Soviet 
attack  and  ballistic  missile  submarines  to  a  more  localized  focus  towards 
conventionally  powered  submarines  operated  by  Third  World  countries. 

Apart  from  the  United  States  and  the  Commonwealth  of  Independent 
States  (CIS),  over  40  other  countries  now  operate  more  than  400  submarines 
worldwide  (Withers,  1992),  half  of  which  are  operated  by  Third  World 
countries.  Some  of  the  more  modern  submarine  types  held  in  Third  World 
inventories  include:  the  Ex-Soviet  KILO-Class,  the  German  Type  209,  the 
British  OBERON-Class  and  the  French  AGOSTA-Class.  Historically,  diesel 
electric  submarine  (SSK)  operations  have  been  limited  in  both  endurance  and 
range.  However,  recent  developments  in  propulsion  technology,  including  air 
independent  propulsion  (AIP),  have  substantially  reduced  snorkeling 
requirements  and  have  also  produced  a  quieter,  more  formidable  SSK.  Due 
to  their  smaller  size,  SSKs  offer  lower  target  strengths  than  their  nuclear 
counterparts;  they  offer  less  of  a  target  to  ping  on  and  thus  produce  less 
return.  When  operating  slowly,  they  are  non-cavitating  and  provide  little 
Doppler  retxmn.  The  net  effect  of  these  improvements  is  to  reduce  the 
detection  opportunities  for  ASW  forces.  In  order  to  gain  control  of  the  seas, 
whether  den5dng  it  to  an  adversary  or  protecting  it  for  self  use,  the  U.  S. 
Navy  must  develop  the  ability  to  conduct  successful  ASW  operations  within 
the  littoral  environments. 
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B.  BACKGROUND 

Shallow  water  acoustics  has  long  been  a  topic  of  ASW  research  and 
has  been  thoroughly  investigated  both  theoretically  and  experimentally  (e.g., 
Akal,  1980;  Frisk  et  al.,  1984;  Rajan  et  al.,  1987).  However,  the  accumulation 
of  theory  and  direct  measurements  has  failed  to  provide  consistent,  accurate, 
and  qualitative  prediction  of  acoustic  propagation  in  shallow  water 
environments.  The  reason  for  this  is  due  almost  entirely  to  the  complexity  of 
the  problem.  In  shallow  water  the  acoustic  medium  consists  of  properties 
that  vary  both  spatially  and  temporally  on  rather  short  scales.  Since  the 
spatial  variation  of  geoacoustic  properties  of  the  bottom/sub-bottom  is  large, 
it  is  unhkely  that  a  complete  and  accurate  data  base  of  bottom/sub-bottom 
geoacoustic  properties  will  ever  be  collected  for  all  shallow  water  areas  of 
interest.  Therefore,  the  U.  S.  Navy  critically  needs  and  accurate,  cost 
effective,  and  practical  method  to  measime  or  deduce  the  geoacoustic 
properties  of  the  sub-bottom.  The  aim  of  this  thesis  is  to  solve  this 
challenging  problem  vital  to  the  U.  S.  Navy’s  abihty  to  accurately  model 
acoustic  propagation  in  shallow  water. 

Unhke  the  deep-water  environment,  acoustic  propagation  in  shallow 
water  is  dominated  by  repeated  interactions  with  the  seafloorand  sub-bottom. 
A  distinguishing  characteristic  of  shallow  water  propagation  is  that  the 
sound-speed  profile  (SSP)  is  generally  nearly  constant  with  depth  or  is 
downward  refi;’acting,  meaning  that  long-range  propagation  takes  place 
exclusively  via  bottom-interacting  paths.  The  deep-water  concept  of  direct 
path  (DP)  and  convergence  zone  (CZ)  acoustic  propagation  does  not  apply  in 
the  shallow  water  environment.  Instead,  the  more  important  propagation 
paths  are  of  the  form  of  refracted  bottom-reflected  or  surface-reflected 
bottom-reflected  propagation.  Generally,  the  high  attenuation  resulting  firom 
repeated  acoustic  interaction  with  the  bottom  and  Hmited  water  depths 
severely  impacts  propagation,  thereby  reducing  detection  ranges  in  shallow 
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water.  Figiire  1  gives  an  example  of  the  degree  of  transmission  loss  (TL) 
variability  as  a  function  of  range  in  shallow  water.  The  plot  illustrates 
average  TL  versus  range  for  shallow  regions  (100-200  m)  around  the  world 
(Urick,  1979).  Two  features  are  of  immediate  interest  in  this  figure.  One  is 
the  spread  of  the  TL  data  in  excess  of  50  dB  at  100  km  due  primarily  to  the 


Figure  1.  Transmission  loss  variability  in  shallow  water  (from  Urick,  1979). 


varying  geoacoustic  makeup  of  the  seafloor.  The  second  feature  is  that 
propagation  is  generally  better  than  that  experienced  with  spherical 
spreading  (20  log  r)  at  short  and  intermediate  ranges  but  is  larger  for 
extended  ranges.  This  feature  is  attributed  to  reduced  spreading  loss  due  to 
the  trapping  of  acoustic  energy  in  the  shallow  water  wave  guide  which  tends 
to  decrease  TL  at  shorter  ranges.  However,  due  to  increased  interaction  with 
the  boundaries,  TL  increases  at  longer  ranges  (Jensen  et  al.,  1994).  The 
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cross-over  range  depends  on  water  depth,  acoustic  jfrequency,  and  the  nature 
of  the  sea  bottom. 

At  low  frequencies  acoustic  propagation  in  shallow  water  is  highly 
dependent  on  the  geoacoustic  properties  of  the  sea  bed.  As  sound  propagates 
through  the  sea  bed,  the  acoustic  pressure  field  is  influenced  by  the  unique 
physical  properties  of  the  sedimentary  sequence  through  which  it  travels  and 
also  by  the  nature  of  the  hthologic  boundaries.  Thus,  an  area  of  considerable 
interest  in  underwater  acoustics  is  the  determination  and  qualitative 
description  of  these  geoacoustic  parameters  that  effect  acoustic  propagation. 
The  traditional  approach  for  describing  the  acoustic  interaction  with  the 
bottom  has  been  to  develop  a  geoacoustic  model  (direct  method)  by  methods 
described  in  Hamilton  (1980).  Provided  adequate  geologic  and  geophysical 
data  exist  in  the  study  area,  this  method  has  proved  invaluable  to  acoustic 
propagation  modeling  efforts. 

A  second  approach  involves  the  inversion  of  these  properties  from  a 
measured  sound  pressure  field.  Numerous  techniques  have  been  cited  in  the 
hterature  for  obtaining  inverse  geoacoustic  solutions.  Frisk  et  al.  (1986)  used 
amplitude  versus  range  data  obtained  using  a  deep-towed  pulsed  continuous 
wave  (CW)  source  and  two  receivers  anchored  near  the  bottom  to  infer 
geoacoustic  properties  of  the  upper  sedimentary  layers  in  deep  water.  The 
inferred  geoacoustic  model  was  derived  from  the  measured  data  using  an 
iterative  forward  modeling  technique.  The  method  rehes  on  time-gating  to 
separate  the  bottom  reflected  arrivals  from  the  surface  reflected  signal.  Zhou 
and  Zhang  (1987)  deduced  geoacoustic  properties  through  normal  mode 
filtering  techniques.  Their  technique  utilizes  dispersion  analysis  and  normal 
mode  measurements  to  determine  acoustic  attenuation  and  propagation 
speeds  for  a  horizontally  stratified  bottom  as  a  function  of  frequency.  This 
type  of  experiment  requires  the  use  of  a  vertical  array  of  hydrophones  to 
successfully  separate  the  modal  arrivals. 
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In  this  thesis  we  use  a  perturbative  inversion  technique  (Lynch  et  al., 
1987;  Rajan,  1992)  to  obtain  the  compressional  sound  speed  as  a  function  of 
depth  for  the  upper  sedimentary  layers  using  as  input  data  the  eigenvalues 
of  the  propagating  sound  field.  The  method  consists  of  measuring  the 
magnitude  and  phase  of  the  pressure  field  as  a  function  of  range  generated 
by  a  CW  point  source  and  numerically  Hankel  transforming  these  data  to 
obtain  the  depth-dependent  Green’s  function  versus  horizontal  wavenumber. 
The  resulting  Green’s  function  contains  information  about  the  discrete  and 
continuous  modal  spectra  of  the  waveguide;  specifically  the  prominent  peak 
in  the  Green’s  function  occxirs  at  wave  numbers  corresponding  to  the 
eigenvalues  for  the  trapped  modes  of  the  waveguide.  Although  the 
sedimentary  layers  have  enough  rigidity  to  support  shear  wave  propagation, 
Fryer  (1978)  indicates  that  the  effect  of  compressional  to  shear  wave 
conversion  at  the  layer  interfaces  is  minimal  for  frequencies  above  20  Hz. 
Thus  in  this  work,  the  sediments  are  treated  as  a  fluid  extension  of  the  water 
column  and  the  shear  wave  effects  are  neglected.  Furthermore,  the  water 
column  is  characterized  as  isospeed  bounded  at  the  sxirface  and  bottom  by  a 
horizontally  stratified  media. 

C.  PROJECT  GEMINI 

A  series  of  narrow-band  bottom  interaction  experiments,  collectively 
called  Project  GEMINI,  were  performed  in  the  Gulf  of  Mexico  near  Corpus 
Christi,  Texas  in  September  1985.  Project  GEMINI  was  designed  to  be  a 
“basehne”  experiment  to  evaluate  low-frequency  acoustic  interaction  with 
the  seabed  (Lynch  et  al.,  1991).  The  intent  of  the  experiment  was  to  obtain 
high-resolution  acoustic  measurements  at  a  relatively  benign  (i.e., 
geologically  simple,  range-independent)  acoustic  environment  which  could 
later  be  used  to  test  or  benchmark  state-of-the-art  acoustic  propagation 
models.  In  order  to  properly  address  the  performance  and  robustness  of 
propagation  models  in  shallow  water,  the  models  must  first  be  evaluated  in 
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regions  where  the  physical  parameters  (e.g.,  sound  speed,  geologic  structure, 
etc.)  affecting  acoustic  propagation  are  well  understood.  A  propagation  model 
that  does  not  perform  well  in  a  geologically  simple  area  is  less  likely  to 
perform  well  in  more  geologically  complex  regions.  The  GEMINI  region  is 
characterized  by  a  smooth,  gently  sloping  seafloor  overl5dng  a  sequence  of 
horizontally-layered  sediments.  Physical  characteristics  of  the  GEMINI 
seafloor  such  as  sediment  thickness,  speed,  and  density  are  well-known  and 
documented  (Matthews  et  al.,  1985).  Thus,  the  experiment  provided  the 
opportunity  to  qualitatively  determine  how  well  one  could  predict  measured 
acoustic  data  using  propagation  models  with  high  quaHty  geoacoustic 
information  as  input  to  the  models. 

Project  Gemini  consisted  of  five  synthetic  aperture  array  experiments 
occup5dng  three  different  array  locations  along  the  continental  shelf  of 
Texas.  The  experiments  consisted  of  towing  a  narrow-band  CW  source  at  a 
fixed  depth  with  output  tonals  of  50  Hz  and  140  Hz,  away  from  a  pair  of 
moored  receivers  (Lynch  et  al.,  1991).  The  track  lines  in  each  of  the 
experiments  extended  from  zero  range  at  the  array  locations  to  a  maximum 
distance  of  about  5000  m.  The  fines  were  oriented  such  that  they  transversed 
upslope  from  the  array  positions  nearly  perpendicular  to  the  isobath  fines. 
The  two  moored  receivers  were  placed  midway  in  the  water  column  and 
approximately  1.5  m  above  the  bottom,  respectively.  A  complete  fisting  of 
array  locations,  source  and  receiver  combinations,  and  dates  of  the 
experiments  are  presented  in  Table  1. 

The  primary  array  site  selected  for  Project  GEMINI,  here  referred  to 
as  the  Rubano  Site,  was  previously  studied  by  the  Applied  Research 
Laboratory,  Pennsylvania  State  University  (ARL/PSU)  in  1975  (Rubano, 
1980).  Two  alternate  sites,  a  shaUow-site  to  the  north  and  a  deep-site  to  the 
southeast  were  also  incorporated  into  the  exercise  (Figure  2).  The  water 
depth  for  the  shallow  site  is  approximately  one-half  that  of  the  Rubano  site 


6 


whereas  the  deep-site  is  about  twice  the  depth.  Prior  to  the  exercise  a 
detailed  environmental  assessment  of  the  GEMINI  region  was  carried  out  by 
researchers  at  the  Naval  Research  Laboratory-Stennis  Space  Center 
Mississippi  (NRL-SSC)  (Matthews  et  al.,  1985).  This  environmental 
assessment  produced  an  acoustical  evaluation  of  the  test  site  which  used  as 
input  a  high-resolution  geoacoustic  model.  The  geoacoustic  model  was 
constructed  from  the  most  current  multi-channel  seismic,  bathymetric,  and 
drill-hole  information  available  at  the  time.  The  geoacoustic  model  was 
specifically  tailored  for  the  Rubano  site  although  it  was  hoped,  given  the 
proximity  of  the  array  locations  and  the  presumed  benign  geologic  properties 
of  the  region,  that  the  model  could  be  extrapolated  to  the  two  other  sites. 


Experiment 

Number 

Date 

Water 

Depth  (m) 

Source  Depth 

(m) 

Receiver 

Depths  (m) 

Description 

1 

9/08/85 

30 

9 

15  and  29 

Rubano  site 

2 

9/09/85 

30 

23 

15  and  29 

Rubano  site 

3 

9/10/85 

62 

46 

32  and  61 

4  1 

21 

9 

16  and  20 

Shallow  Site 

5 

9/12/85 

62 

9 

32  and  61 

llllllllllEfH^ 

Table  1.  Description  of  Project  GEMINI. 
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Figure  2.  Geographic  locations  of  the  Project  GEMINI  experiments  off  the  Texas 
Coast.  Isobaths  in  meters. 

D.  OBJECTIVE 

In  shallow  water  the  geologic  and  geoacoustic  properties  of  the  bottom 
can  have  significant  vsiriations  over  rather  small  ranges.  Current 
geoacoustic  modeling  methods  (point  models)  only  address  the  physical 
properties  of  the  sediment  and  rock  as  a  function  of  depth  below  the  sea  floor 
at  a  single  point  and  not  as  a  function  of  horizontal  distance.  In  order  to 
adequately  accoxint  for  the  spatial  variations  in  geologic  characteristics  a 
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series  of  geoacoustic  models  must  be  developed  to  account  for  this  variability. 
However,  as  stated  previously,  such  spatially  varying  data  are  not  generally 
available.  An  alternative  approach  is  to  infer  these  geoacoustic  parameters 
through  geoacoustic  inversion  techniques. 

The  objectives  of  this  study  are  three-fold.  The  first  objective  is  to 
compare  several  propagation  model  estimates  of  TL  with  measured  low- 
frequency  acoustic  propagation  data  collected  during  Project  GEMINI.  Such 
a  comparison  will  determine  which  of  several  TL  models  most  closely 
approximates  the  observed  TL  pattern,  i.e.,  which  model  solves  the  forward 
problem  best.  The  acoustic  models  will  have  as  input  location- specific,  high- 
resolution  geoacoustic  characterizations  of  the  seafloor.  The  acoustic  models 
that  wiU  be  used  in  this  study  are  the  SACLANTCEN  normal  mode  model 
(SNAP),  the  NRL  normal  mode  model  (KRAKEN),  and  the  NRL  finite 
element  parabolic  equation  model  (FEPE). 

The  second  objective  is  to  infer  the  geoacoustic  properties  of  the 
seafloor  in  the  GEMINI  region  using  the  perturbative  inversion  method.  The 
underlying  assumption  of  perturbative  inversion  theory  is  that  a  reasonable 
“first  guess”  of  the  unknown  geoacoustic  properties  must  be  incorporated  as  a 
starting  point  in  the  computation  process  (Frisk  et  al.,  1989).  Using  a 
detailed,  location-specific  geoacoustic  model  as  a  starting  point  in  the 
inversion  process,  it  will  be  shown  that  a  geoacoustic  model,  derived  through 
inversion,  results  in  improved  TL  model  estimates  over  those  presented 
earlier. 

Lastly,  an  objective  vitally  important  to  the  Navy’s  operational 
performance  in  shallow  water  will  be  pursued.  In  many  shallow  water 
regions  little  geologic  or  geophysical  information  is  readily  available  to 
formulate  useful  geoacoustic  characterizations  of  the  bottom.  If,  however  an 
“educated  first  guess”  for  a  geoacoustic  model  is  made  in  a  shallow  water 
area  of  interest,  it  is  important  to  determine  how  accurate  the  perturbative 
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inversion  technique  will  be  in  estimating  the  actual  geoacoustic  properties. 
Moreover,  the  TL  model  results  using  this  process  must  be  assessed  to  see 
how  comparable  they  are  to  those  formulated  using  a  high  quahty 
geoacoustic  model  as  an  input. 
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II.  NORMAL  MODE  SOLUTIONS  IN  A  WAVEGUIDE 


A.  MODAL  PROPAGATION 


The  inverse  method  considered  in  this  thesis  requires  as  input  data 
the  measured  horizontal  wave  numbers  or  eigenvalues  which  characterize 
the  modal  propagation  of  sound  through  the  ocean  and  seafloor.  Since  the 
ideas  presented  in  normal  mode  theory  are  essential  to  the  inversion  process, 
a  brief  discussion  of  normal  mode  propagation  is  presented.  Normal  mode 
theory  provides  a  ‘full  wave’  solution  of  the  wave  equation  in  the  ocean 
waveguide  which  takes  into  account  the  physical  characteristics  of  the 
medium,  frequency  dispersion  effects,  and  interaction  of  the  sound  field  at 
the  boundaries  (Casey,  1988). 

The  propagation  of  monochromatic  sound  in  the  ocean  is  described 
mathematically  by  solutions  to  the  reduced  wave  or  Helmholtz  equation. 
Using  the  assumption  of  cylindrical  symmetry,  the  spatial  part  of  the 
acoustic  pressure  p(r,z,Zo),  due  to  a  point  source  at  range  r  =  0  and  depth  z  = 
Zo  with  a  sinusoidal  or  harmonic  time  dependence  exp(-i(Dt),  satisfies  the 
Helmholtz  equation 


1 

r 


k^z) 


p(r,z,  zj 


(1) 


where  p(r,z,Zo)  represents  the  acoustic  pressure  at  range  r  and  depth  z,  <o  is 
the  radian  frequency,  and  k(z)  =  G)/c(z)  is  the  total  wave  nrunber. 

Appropriate  boundary  conditions  must  be  specified  in  order  to 
determine  a  unique  solution  to  Equation  (1).  At  the  interfaces  where  there 
are  significant  variations  in  properties,  e.g.,  at  the  water  sediment  interface, 
continuity  of  pressure  and  normal  particle  velocity  must  be  met.  The 
symmetry  and  range  independence  of  Equation  (1)  leads  to  the  method  of 
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separation  of  variables  in  solving  the  equation.  By  expressing  the  pressure 
field  p(r;z,z,>)  by 

p(r,z,Zo)  ==Z(z)R(r)  (2) 

and  substituting  this  into  Equation  (1)  leads  to  the  following  expression 


The  first  term  in  Equation  (3)  is  a  function  of  r  only  whereas  the 
second  term  is  a  function  of  z  only.  This  equality  is  satisfied  if  and  only  if 
both  expressions  are  equal  to  the  same  constant  of  proportionality  which 
leads  to  the  two  separated  equations: 


d^R(r)  1  dR(r) 

2  +  j  +k,R(r)-0 

dr  r  dr 

(4) 

and 

dz 

(5) 

where  (kr)^  and  (kz)^ 

are  the  separation  constants  defined  by 

k^  =k"(z)-k^ 

(6) 

Equation  (4)  is  a  function  of  range  only,  and  thus  from  Equation  (6)  k,  ,  the 
horizontal  component  of  the  wave  number,  must  be  constant  throughout  the 
waveguide.  This  can  also  be  expressed  in  terms  of  Snell’s  law  since 

G) 

=  ksin6  =  — sin6  (7) 

c 

where  9  is  the  vertical  angle  of  incidence  and  to  is  constant,  and  (sin0)/c  is 
constant  (Clay  and  Medwin,  1977).  Accordingly,  at  a  zero  incidence  angle 
(0  =  0°),  there  is  no  horizontal  contribution  to  the  total  wavenumber.  The 
solution  to  Equation  (4)  is  a  cylindrical  Bessel  or  Hankel  function  of  the 
zeroth  order. 

R(r)  =  Hi^>(k,r)  (8) 

For  distances  of  more  than  a  few  wavelengths  from  the  source,  the 
asymptotic  expansion  of  H|,‘^(k,r)  (Haberman,  1987)  can  be  expressed  as 
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(9) 


'  Ttk^r 


C09 


7T^ 

4) 


The  exponential  expansion  of  the  cosine  term  in  Equation  (9)  yields  the 


following  characteristic  equation 


H“’(k,r) 


^  (  ■  \ 
•  1  Irt 

exp  ik  r-—  +exp  -ik^r  +  — 

^  4;  ^  4; 


(10) 


When  associated  with  a  harmonic  time  dependent  term,  Equation  (10) 
represents  the  total  wave  motion  consisting  of  an  inward  propagating  wave 
represented  by  the  first  term  on  the  right  side  of  Equation  (10)  and  an 
outward  propagating  wave  represented  by  the  second  term  on  the  right  of 
Equation  (10).  Since  we  are  only  interested  in  an  outward  propagating, 
exponentially  decaying  waveform,  only  the  second  term  is  used  in  the 
formulation  (Clay  and  Medwin,  1977). 


Hl^k.r) 


(11) 


1.  Rigid  Bottom  Solution 

Assuming  that  the  ocean  surface  represents  a  pressure  release 
surface  and  that  the  ocean  bottom  at  z  =  h  is  perfectly  rigid,  allowing  no 
acoustic  penetration,  the  “normal  mode”  differential  equation  (Equation  (5)) 
and  the  accompanying  boundary  conditions  are: 

^^  +  k;z<z)=0  (12) 

Z(0)=0,  ^^  =  0  (13) 

By  solving  this  equation  as  a  Sturm-Liouville  eigenvalue  problem,  the 
general  solution  to  Equation  (12)  is  of  the  form 

Z(z)=Cisin(kzz)+C2COs(kzz)  (14) 

where  Ci  and  C2  are  arbitrary  constants  and  ks  is  the  vertical  wavenumber. 
The  boundary  condition  at  the  free  surface  (z=0)  requires  that  the  acoustic 
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pressure  to  be  zero  for  all  range,  r,  and  time,  t,  thus  requiring  C2=0. 
Furthermore,  an  incident  plane  wave  on  the  sea  surface  is  assumed  to  be 
completely  reflected,  undergoing  a  180°  phase  shift  upon  interaction  with  the 
free  surface.  The  rigid  bottom  boundary  condition  is  characterized  by  the 
constraint  that  the  acoustic  pressure  is  a  maxima  at  the  seafloor  at  z=h  and 
that  the  normal  particle  velocity  vanishes  on  the  surface  resulting  in 
complete  reflection  (Casey,  1988).  This  means  that  kzZ  =  it/2,  3:t/2,  57t/2...  or 
more  generally 

(kz)m  =  (3t/h)(2m-l),  m  =  l,  2,  3...  (15) 

(kz)m  is  termed  an  eigenvalue  and  the  solution  is  an  eigenfunction  (normal 
mode)  of  the  problem.  There  exists  an  infinite  number  of  eigenvalues  and 
thereby  an  infinite  number  of  eigenfunctions  Zm(z)  as  solutions  to  this 
boundary  value  problem.  Each  solution  is  characterized  by  a  mode  shape 
Zni(z)  and  a  vertical  propagation  constant  (kz)m.  The  subscript  m  designates 
the  modal  number  for  a  particular  eigenvalue  or  eigenfunction. 

The  pressure  field  in  Equation  (2)  can  now  be  written  as 

p(r,z)  =  £R,„(r)Z^(z)  (16) 

m=l 


Substituting  this  expression  for  pressure  into  Equation  (1)  yields: 

-5(r)5(z-z„) 


a 


‘l  d 

_r  dr 

1  dr  j 

Z.+k:R.Z 


z  m 


27tr 


(17) 


The  eigenfunctions  of  the  Sturm-Liouville  problem  form  an  orthonormal  set 
of  square  integrable  functions  on  0<h<z  such  that 


n 

I 


_jOm*n 


p(z) 


1  m  =  n 


(18) 


Multiplying  both  sides  of  Equation  (17)  by 

K(z) 


P(z) 


■dz 


(19) 


and  applying  the  orthogonality  relationship  of  (18)  yields  the  following 
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Idf  dR„ 


r  drv  dr 


+  k?  R  = 


-8(r)Z„,(zJ 

27trp(zJ 


The  solution  to  Equation  (20)  is  a  Hankel  function  of  the  first  kind  such  that 

Rm(r)  =  ^^^Z,„(zJH<'Hk,.,r)  (21) 

Substituting  the  expression  for  H™(k^r)(knr)  found  in  Equation  (11),  the 
solution  to  Equation  (1)  can  be  expressed  as  a  sum  of  normal  modes 

i  ^  exp(ik,  r) 

p(r,z,zj=  exp(^)^  Z„(z)Z„,(zJ - j=^  (22) 

p(zjv87ir 


Equation  (22)  describes  the  pressure  field  as  a  linear  superposition  of 
traveling  cylindrical  waves  propagating  outward  from  a  source.  Each  of 
these  modes  is  a  standing  wave  in  the  vertical  direction  with  a  unique  depth 
distribution  described  by  the  eigenfunction  Zm(z)  and  the  vertical  modal 
wavenumber  (kz)m  The  vertical  wavenumber  must  ensure  that  a  pressure 
node  exists  at  the  surface  and,  in  the  case  of  a  perfectly  rigid  bottom,  a 
pressure  anti-node  at  the  seafloor.  The  eigenfunctions  Zm(z)  essentially 
determine  the  depth  dependence  of  the  acoustic  pressure  field.  For  a  given 
water  depth  and  source  frequency,  only  a  finite  number  of  modes  may  exist. 
For  large  mode  numbers  rn,  the  horizontal  wavenumber  becomes  imaginary 
and  the  solution  defined  by  Equation  (22)  rapidly  decays  with  increased 
horizontal  ranges. 


2.  Penetrable  Bottom  Solution 

If  the  ocean  bottom  is  not  a  perfect  reflector,  sound  energy  will 
penetrate  into  the  sea  bed  and  solutions  to  the  Helmholtz  equation  must  be 
found  both  in  the  water  column  and  the  ocean  bottom.  The  solution  to 
Equation  (1)  for  the  penetrable  bottom  case  can  be  obtained  by  utilizing  the 
Green’s  function  technique.  The  conditions  imposed  upon  this  case  are  the 
same  as  that  of  the  rigid  bottom  with  the  exception  that  the  velocity 
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potential  is  discontinuous  at  the  seafloor.  A  brief  derivation  of  the  Green’s 
function  is  now  given. 

Assuming  a  constant  density  water  column  of  thickness  (h),  sound 
speed  c(z),  bounded  at  the  surface  and  bottom  by  a  horizontally  stratified 
medium  (Figure  3),  the  spatial  part  of  the  sound  pressure  field  (p)  due  to  a 
harmonic  time  dependent  point  source  located  at  r  =  0  and  z  =  Zo  satisfies 
Equation  (1). 

The  solution  to  Equation  (1)  can  now  be  written  as  the  zero-order 
Hankel  transform  of  the  depth-dependent  Green’s  function  g(kr,z).  This  is 
achieved  by  taking  the  zero-order  Hankel  transform  of  both  sides  of  Equation 
(1)  where  the  Hankel  transform  pair  is  defined  by  the  following: 


p(r,z,)  =  Jg(k,,z)J„(k,r)k,dk, 

0 

(23) 

CO 

g(k,,z)  =  Jp(r)J„(k^r)rdr 

(24) 

0 


The  Hankel  transform  variable,  kr,  transforms  the  solution  from  a 
depth  and  range  space  p(z,r)  to  a  solution  as  a  function  of  depth  and 
horizontal  wave  number  space  (z,  kr)  (Evans,  1975). 


Surface  Plane  Wave  Reflection  Coefficient  R  ^ 

i  ^ 

z 

•  Source 
(O.Zol 

•  Receiver 
(r,  z) 

_ Xcb _ 

Bottom  Plane  Wave  Reflection  Coefficient  R  ^ 

Horizon taiiv  Stratified  Bottom 

Figure  3.  Horizontally  stratified  ocean  model. 
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After  applying  Equation  (24)  to  Equation  (1)  and  re-ordering  terms, 
Equation  (1)  is  transformed  into  the  depth  separated  wave  equation 

^+k^(z)-k,  g[k  2,z  )=-25(z-zJ  (25) 

Here  g(kr,z)  is  the  depth-dependent  Green’s  function,  Jo  is  the  zero-order 
Bessel  function,  r  is  horizontal  range,  and  kr  is  the  horizontal  wave  number. 
The  general  solution  to  Equation  (25)  can  he  written  as  the  sum  of  the 
particular  solution  (non-homogeneous  part)  and  the  characteristic  solution 
(homogenous  part).  The  characteristic  solution  of  Equation  (25)  is  of  the  form 

g(kr)  =  (t»S+(t)B  (26) 

where  and  (|)b  are  linearly  independent  solutions  of  Equation  (25)  chosen 
such  that  bs  satisfies  the  surface  boundary  condition,  and  <t)B  satisfies  the 
bottom  boundary  condition.  By  applying  the  method  of  variation  of 
parameters,  the  Green’s  function  can  be  defined  by: 

O.z.z.  (27) 

W(zJ 


0  <  z  <  z„ 


g(kr,z)  = 


-2(t>s(l<^r.Zo)<l)B(k,,Z) 

W(zJ 


z„  <  z  <  h 


such  that  W(zo)  is  the  Wronskian  given  by 

W(Zo)  =  <|)b(Zo)  <1)'s(Zo)-  (I)'b(Zo)  (t)s(Zo)  (29) 

where  the  prime  indicates  differentiation  with  respect  to  z.  (|)s  and  (j)B  also 
must  satisfy  the  following  impedance  relations. 

atz=  0 


atz  =  h 


^  and  ^8  are  input  quantities  that  contain  boundary  conditions  necessary  for 
the  calculation  of  the  Green’s  function  (Frisk  et  al.,  1980).  Both  quantities 
are  dependent  on  the  horizontal  wavenumber,  frequency,  and  the  acoustic 
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properties  of  the  bottom  and  the  surface,  depends  on  the  water  density, 
Pvv,  and  the  sound  speed  c(z),  at  the  water-bottom  interface,  while  depends 
on  fHv  and  the  sound  speed  c(z)  at  the  water-surface  interface.  Equation  (30) 
therefore  can  be  expressed  in  terms  of  plane  wave  reflection  coefficients,  Rb 
and  Rs.  for  the  half-space  corresponding  to  sound  speed  cb  and  cs  bounding 
the  bottom  and  the  surface,  respectively. 

(1)b(z)  =  Ci[exp(-ikzz)+  RB(kr)  exp(ikzz)]  (31) 

(t)s(z)  =  C2[exp(-ikzz)-i-  R8(kr)  exp(ikzz)]  (32) 

where  kz^+kr^  =k2=(a)/c)2  is  the  total  wavenumber  and  Ci,  C2  being  arbitrary 
constants.  Equation  (31)  can  also  be  re-written  as 


= 


*1^B 


RI  +  Rb) 

kRl-Re)’ 


at  z  =  0 


-i  l  +  Roexp(-2ik7h)  , 

- 2 ^  ,  at  z  =  h 

kz  l_l-Rsexp(-2ik2h) 


(33) 

(34) 


Substituting  these  expressions  into  Equations  (27)  and  (28)  and 
further  assuming  that  the  ocean  surface  is  a  pressure  release  boundary  (i.e., 
Rs  =  -1),  the  following  expression  for  the  depth  dependent  Green’s  function 


can  be  found: 

g(kr,z)  = 


(35) 


where  z+  =  z-t-Zo,  z-  =  z-Zo 


In  order  to  evaluate  the  integral  expression  for  p(r,z)  the  Bessel 
function,  Jo(k,r)  in  Equation  (23),  is  expressed  in  terms  of  Hankel  functions 

as; 
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■lo(k,r)  =  5[Hi;'(k,r)  +  H®(k,r)] 


(36) 


The  following  relation 

exp(-i3Tk,r)  =  -H[,'Hk,r)]  (37) 

allows  the  range  of  integration  of  the  transform  integral  to  be  evaluated  over 
the  entire  real  kr  axis  from  -(»  to  oo  (Casey,  1988).  The  solution  of  Equation 
(23)  is  then  found  through  contour  integration  methods  (Evans,  1975)  and  is 
expressed  as  the  sum  of  discrete  and  continuous  portions  of  the  normal  mode 
spectrum  given  by 


p(r,z)  =  m 


? 


ZJz)Z„(zJH<'>(k,j)  +  I(r) 


where  the  eigenvalues  km  and  eigenfunctions  Zm  satisfy 


+  k'(z)-k' 
V 


Z.(z)  =  0 


(38) 


(39) 


The  discrete  sum  in  Equation  (38)  corresponds  to  modes  perfectly 
trapped  in  the  waveguide  and  is  of  the  same  form  as  that  of  the  rigid  bottom 
solution.  Typically,  this  trapped  mode  sum  dominates  the  long  range 
behavior  of  the  acoustic  field  whereas  the  continuum  contribution  I(r)  is  only 
significant  at  short  ranges  (Lynch  et  al.,  1987).  The  continuum  is  rapidly 
attenuated  with  increasing  range  fi’om  the  source  and  is  typically  neglected 
for  ranges  over  a  few  wavelengths  fi:om  the  source.  The  eigenvalues  for  the 
penetrable  bottom  case  have  the  same  interpretation  as  those  in  the  rigid 
bottom  solution.  However,  they  will  take  on  different  values  since  sound  is 
propagating  in  both  the  water  column  and  the  sub-bottom  and  the 
characteristic  equation  will  be  different. 

The  modal  representation  is  closely  correlated  with  the  analytic 
properties  of  the  Green’s  function.  A  knowledge  of  the  behavior  of  g(kr)  (for 
example,  the  number  of  resonances  and  their  spatial  location  in  regard  to 
horizontal  wavenumber)  provides  information  about  the  nature  of  the  modal 
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propagation  in  the  waveguide.  The  Green’s  function  is  characterized  by  a 
finite  valued  continuum  with  a  discrete  set  of  resonances  occurring  when  the 
denominator  of  Equation  (34)  is  set  equal  to  zero.  These  poles  of  g(kr,z) 
correspond  to  the  eigenvalues  of  the  trapped  and  virtual  modes  excited  in  the 
waveguide.  The  position  and  magnitudes  are  highly  influenced  by  the 
inherent  acoustic  properties  of  the  bottom  (Frisk  et  al.,  1986)  and  the 
geometry  of  the  waveguide.  As  a  result,  the  basic  principle  underlying  the 
inverse  method  presented  in  this  paper  is  first  to  obtain  an  estimate  of  g(kr) 
from  measurements  of  p(r)  by  numerically  performing  the  Hankel  transform 
in  Equation  (23).  The  modal  features  of  g(kr)  are  then  used  to  determine  a 
geoacoustic  model.  A  specific  perturbation  technique  relating  the 
eigenvalues  of  the  trapped  modes  to  the  geoacoustic  model  are  given  in  the 
next  section. 

B.  PERTURBATIVE  INVERSION  TECHNIQUE 

A  detailed  discussion  of  the  perturbative  inversion  theory  that  is  the 
basis  for  the  analysis  in  this  thesis  is  presented  in  Rajan  et  al.  (1987).  The 
objective  of  the  perturbative  inversion  technique  is  to  estimate  the  variation 
in  geoacoustic  parameters  (i.e.,  sound  speed,  density,  attenuation,  etc.)  from 
the  differences  between  the  measured  eigenvalues  and  the  eigenvalues 
computed  for  a  background  geoacoustic  model  (Frisk  et  al.,  1986).  In  the 
perturbative  technique,  an  initial  background  geoacoustic  model  is  assumed 
and  the  differences  of  the  geoacoustic  parameters  from  the  geoacoustic 
parameters  in  the  initial  background  model  are  calculated.  For  the  purpose 
of  this  thesis,  we  are  concerned  with  determining  only  the  compressional 
sound  speed  profile  in  the  bottom,  c(z),  for  the  deep  and  shallow  GEMINI 
sites.  Perturbative  techniques  to  determine  other  geoacoustic  parameters 
(i.e.,  density,  attenuation,  etc.)  are  not  addressed.  However,  these 
parameters  can  be  modified  to  be  geologically  consistent  with  the  final 
estimate  of  c(z).  As  a  first  step,  an  initial  background  geoacoustic  model  for 
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c(z)  is  assumed,  and  a  major  objective  of  this  thesis  is  to  evaluate  the 
perturbative  technique  when  the  initial  estimate  of  c{z)  must  be  based  on 
little  to  no  geoacoustic  information  about  the  area.  Perturbation  theory  then 
provides  the  means  for  determining  the  difference  Ac(z)  between  the 
background  model  and  the  true  earth  model. 

From  first-order  perturbation  theory,  an  integral  equation  relating  the 
perturbation  of  eigenvalues  is  given  by: 


T  Z  ^k(z)^ 

where  Akm  represents  the  difference  between  the  experimentally  measured 
modal  eigenvalues  and  the  eigenvalues  computed  from  the  background 
geoacoustic  model,  km  and  Zm(z)  are  the  mth  eigenvalues  and  normalized 
eigenfunctions  for  some  assumed  background  geoacoustic  model,  p(z)  is  the 
density  profile  for  the  background  geoacoustic  model,  k(z)  is  the  acoustic 
wavenumber  defined  by  co/c(z)  (where  co  is  the  angular  frequency  and  c(z)  is 
sound  speed  profile  of  the  geoacoustic  model),  and  Ac(z)  is  the  difference  in 
sound  speed  between  the  background  and  true  models.  All  terms  in 
Equation  (40)  are  known  except  for  Ac(z)  which  is  the  solution  being  sought. 

Equation  (40)  is  a  Fredholm  integral  equation  of  the  first  kind  for  Ac(z) 
and  generally  does  not  produce  a  unique  solution  unless  constraints  are 
made  upon  its  solution  (Lynch  et  al.,  1987).  Reordering  the  known  quantities 
in  Equation  (40 )  such  that 


results  in 


K(m,z) 


Z.rk(z)^ 

k„PoC„(z) 


(41) 


Ak, 


=  lK(m. 


z)  Ac(z)dz 


(42) 
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K(m,z)  is  defined  as  the  kernel  of  the  integral  equation  and  consists  of  the 
known  physical  parameters  for  a  discrete  mode  and  depth  step.  Using 
methods  proposed  by  Backus  and  Gilbert  (1970),  a  solution  for  Ac(z)  is  given 
by: 

Ac(z)  =  ^^K(m,z)Ak„A''(m,n)  (43) 

m  n 

where  A(m,ri)  is  the  spectral  decomposition  of  the  kernel  K(m,z)  and 
A(m ,  n)  =  I  K(m ,  z)  K  ^  (n,  z)  dz  (44) 

a 

The  solution  given  by  Equation  (43)  is  called  the  Backus-Gilbert 
solution  for  perfect  data.  This  solution  represents  the  most  “Dirac  delta  like” 
resolution  kernel  and  was  used  in  this  analysis  because  smoothing  the 
solution  was  not  desirable  until  after  c(z)  was  estimated.  This  represents  a 
solution  without  constraints  and  is  considered  a  fundamental  property  of  the 
perturbative  technique  not  addressed  originally  in  Raj  an  et  al.,  (1987).  If 
constraints  are  used  in  the  solution  for  c(z),  the  most  natural  (i.e., 
unconstrained)  solution  characteristics  will  not  be  realized.  The  “price  to 
pay”  for  the  natural  or  “unconstrained”  approach  is  instability  in  the  solution 
when  the  ratio  of  the  highest  to  lowest  eigenvalue  is  large.  In  that  instance, 
it  is  well  known  that  the  pseudo-inverse  to  the  matrix  in  Equation  (44)  is  ill 
conditioned,  but  a  stabilization  technique  based  on  the  Frobenious  norm  was 
used.  Also,  a  simple  “box  car”  average  of  c(z)  in  depth  was  performed  to 
dampen  fluctuations  in  c(z)  with  depth. 

As  a  first  step  to  the  inversion  process,  measured  complex  pressure 
field  versus  range  data  was  transformed  from  the  range  domain  to  the 
horizontal  wave  number  domain  (depth-dependent  Green’s  function)  via  the 
Hankel  transform  pair.  As  discussed  earfier,  the  normal  modes  of  the 
waveguide  show  up  as  strong  resonance  peaks.  The  measured  modal 
eigenvalues  are  experimentally  measured  by  simply  picking  the 
peaklocations  in  \g(kr)  t  (Lynch  et  al.,  1991).  Figure  4  is  a  typical  Green’s 
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function  versus  horizontal  wavenumber  plot  for  a  particular  environment.  In 
this  instance,  each  of  the  peaks  (modes  1-4)  has  a  characteristic  wavenumber 
or  eigenvalue  identified  in  our  formulation  as  kmeas. 


Green’s  Function  versus  Horizontal  Wavenumber 
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Figure  4.  Plot  of  the  depth-dependent  Green’s  function  |  g(kr)  |  versus  horizontal 
wavenumher  k  at  140  Hz.  The  normal  modes  of  the  particular  waveguide  are 
identified  by  poles  of  g(kr). 


The  NRL  normal  mode  model,  KRAKEN,  was  used  to  calculate  the 
modal  eigenvalues  and  eigenfunctions,  km  and  Zm(z),  respectively,  using  the 
assumed  background  geoacoustic  model.  In  this  section  we  assume  a 
background  geoacoustic  model  with  a  linear  sound  speed  gradient  and  with 
constant  density  and  attenuation  with  depth  to  illustrate  the  basic  concepts 
of  the  inversion  technique.  In  order  to  demonstrate  the  effectiveness  of  the 
perturbative  inversion  technique,  the  initial  sound  speed  at  the  water-bottom 
interface  and  the  sound  speed  gradient  were  also  varied.  Once  km  and  Zm(z) 
are  calculated,  Akm  is  found  by  the  relation 

Akm  =  kmeas  -  km  (m=l,2,3...)  (45) 


These  values  were  entered  into  Equation  (41)  and  the  result,  Ac(z),  is 
calculated  by  Equation  (43)  and  represents  the  perturbation  of  the  sound 
speed  profile  about  that  proposed  in  the  initial  background  model.  The  sound 
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speed  profile  is  then  adjusted  by  subtracting  Ac(z)  from  the  initial  value  of 
c(z).  The  perturbative  process  continues  until  Ac(z)  is  minimized. 

1.  Synthetic  Test  Cases 

In  this  section  we  will  apply  the  perturbative  inversion  technique  to 
a  set  of  synthetic  data  for  which  the  correct  solution  is  known.  The  synthetic 
earth  model  along  with  the  two  initial  background  geoacoustic  models  are 
shown  in  Figirre  5.  The  synthetic  earth  model  consists  of  a  50  m  thick 
isospeed  water  column  (1500  m/s),  overlying  a  10  m  thick  isospeed  sediment 
layer  (1800  m/s),  which  comformably  hes  over  a  60  m  thick  linear  gradient 
(Vc  =1.50  s-i)  sedimentary  layer.  The  basement  consists  of  an  isospeed  half- 
space  with  a  compressional  wave  speed  of  1900  m/s.  The  density  is  assumed 
to  be  constant  p  =  1.6  g/cm'^  for  aU  sediment  layers  and  half-spaces.  Although 
the  water  column  sound  speed  profile  for  this  example  is  assumed  to  be 
isospeed,  it  can  consist  of  any  structure.  In  general,  the  water  column  profile 
wiU  be  known  fi-om  available  XBT  and  CTD  measurements  and  will  thus 
drop  out  of  the  perturbation  integrals,  leaving  only  the  sediment  properties  to 
be  determined  (Rajan  et  al.,  1987). 


Figure  5.  Properties  of  the  synthetic  earth  and  two  hackgound  model. 
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The  frequencies  that  will  be  considered  in  the  synthetic  examples  are 
25,  50  and  150  Hz.  Sub-bottom  penetration  by  acoustic  energy  is  highly 
dependent  on  frequency  with  lower  frequencies  being  able  to  penetrate 
deeper  into  the  seafloor  than  higher  frequencies.  As  a  result,  lower 
frequencies  are  able  to  resolve  the  deeper  sedimentary  structure,  whereas 
the  higher  frequencies  resolve  the  finer  near-surface  structure.  The  synthetic 
Green’s  function  versus  horizontal  wavenumber  plots  for  each  of  the 
frequencies  are  shown  in  Figures  6a-c  These  poles  of  g(kr,z),  as  discussed 
previously,  correspond  to  the  eigenvalues  of  the  trapped  and  virtual  modes 
that  are  excited  for  the  particular  waveguide  at  the  given  frequency.  In 
general,  higher  frequencies  excite  more  modes  and  thus  have  more  trapped 
eigenvalues  than  do  lower  frequencies.  The  modal  eigenvalues  for  the 
synthetic  earth  model  (kmeas)  are  hsted  in  Table  2. 

Two  initial  geoacoustic  background  models  were  assumed  in  this 
example:  (1)  an  1850  m/s  isospeed  model,  and  (2)  a  Hnear  gradient  model  c(z) 
=  1750  +  2.0  s  Hz)-  For  the  purposes  of  these  test  cases,  density  was  taken  to 
be  a  constant  p  =  1.6  g/cm^,  and  compressional  attenuation  was  assumed 
neghgible.  The  modeled  eigenvalues  (km)  for  the  initial  background  models 
are  shown  in  Tables  3  and  4. 
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Green’s  Function  versus  Horizontal  Wavenumber 
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Figure  6.  Synthetic  Green’s  function  magnitude  versus  wavenumber  for  25  Hz. 


Figure  7.  Synthetic  Green’s  function  magnitude  versus  wavenumber  for  50  Hz. 
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Figure  8.  Synthetic  Green’s  fiinction  magnitude  versus  wavenumber  for  150  Hz. 


Model  Frequency  Mode  Eigenvalue  (kmeas) 


150 _ 3 _ 0.617630422 

150  4  0.583889604 


Table  2.  Modal  eigenvalues  for  the  synthetic  earth  model  (kmeas! 
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Model 

Frequency 

(Hz) 

Eigenvalue 
Iteration  1 

Eigenvalue 
Iteration  2 

Eigenvalue 
Iteration  3 

Eigenvalue 
Iteration  4 

Background 

25 

0.095567562 

0.095840216 

Model  B2 

50 

0.203297004 

0.203447193 

0.203428298 

50 

0.182848796 

0.183509603 

0.183428794 

150 

0.625649512 

0.625688612 

0.625676215 

0.625808477 

150 

0.617533207 

0.617683887 

0.617635310 

0.618125618 

150 

0.603655517 

0.603984475 

0.603875697 

0.604930580 

150 

0.583521008 

0.584122181 

0.583913684 

0.587231219 

Table  3.  Modal  eigenvalues  for  the  initial  background  geoacoustic  model  B1  (km). 


Model 

Frequency 

(Hz) 

Eigenvalue 
Iteration  1 

Eigenvalue 
Iteration  2 

Eigenvalue 
Iteration  3 

Eigenvalue 
Iteration  4 

Background 

25 

0.096081510 

0.095852643 

0.095839143 

Model  B2 

50 

0.203560099 

0.203439906 

0.203429595 

50 

0.183980197 

0.183473706 

0.183432102 

150 

0.625703216 

0.625677884 

0.625686586 

0.617737293 

0.617641628 

0.617674530 

0.604090273 

0.603888690 

0.603957891 

150 

0.584284186 

0.583935380 

0.584054530 

Table  4.  Mod 

[al  eigenvalues  for  the  initial  background  geoacoustic  model  B2  (km). 

Tables  3  and  4  list  the  computed  eigenvalues  for  the  initial 
background  models  for  each  iteration  in  the  inversion  process.  These 
eigenvalues  (km),  as  discussed  previously,  are  subtracted  from  the  measured 
eigenvalues  (Table  2)  to  form  the  data  vector  Ak(m).  Generally,  the  lower 
frequencies  (i.e.,  fewer  modes)  required  fewer  iterations  in  order  for  solutions 
to  converge.  Figures  9-11  display  the  residts  of  the  pertubative  inversion 
technique  to  recover  the  synthetic  earth  sound  speed  profile  from  the  initial 
background  model  Bl. 
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Figure  9.  25  Hz  Inversion  of  synthetic  earth  model  from  initial  background  model 
Bl. 


Figure  10.  50  Hz  Inversion  of  synthetic  earth  model  from  initial  background  model 
Bl. 
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Figure  11.  150  Hz  Inversion  of  synthetic  earth  model  from  initial  background 
model  Bl. 


As  evident  from  the  plots,  the  sound  speed  profiles  are  corrected 
in  the  proper  direction  relative  to  the  background  model.  The  profile  in  the 
upper  tens  of  meters  are  lower  in  speed  relative  to  the  background,  whereas 
the  profile  beyond  about  90  m  depth  is  virtually  unaffected  by  the 
perturbation.  Clearly,  the  25  and  50  Hz  signals  correct  the  profile  to  greater 
depths  than  the  140  Hz  signal.  However,  the  140  Hz  signal  approximates  the 
sound  speed  at  the  water-bottom  interface  more  accurately  than  the  lower 
frequencies.  Also  apparent  from  the  plots  is  that  as  the  normal  mode 
function  decreases  with  depth,  the  inferred  sound  speed  profile  returns  to  the 
background  profile. 

Figures  12-14  show  the  results  of  recovering  the  synthetic  earth 
model  using  the  background  model  B2.  The  results  are  qualitatively  the 
same  as  the  previous  background  model  Bl.  Even  though  the  background 
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models  differ  significantly,  the  solutions  converge  about  a  sound  speed  profile 
which  is  strikingly  similar  to  the  synthetic  earth  model. 


Sound  Speed  (m/s) 

1650  1700  1750  1800  1850  1900  1950 


Figure  12.  25  Hz  Inversion  of  synthetic  earth  model  from  background  model  B2. 


Sound  Speed  (m/s) 

1650  1700  1750  1800  1850  1900  1950 


Figure  13.  50  Hz  Inversion  of  synthetic  earth  model  from  background  model  B2. 
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Figure  14.  150  Hz  Inversion  of  synthetic  earth  model  from  background  model  B2. 


The  source  angles  in  the  data  and  the  models  are  predominantly 
small  (<  20°)  and  thus  the  perturbation  technique  alone  can  not  be  expected 
to  resolve  geoacoustic  features  much  beyond  an  acoustic  wavelength  below 
the  water-bottom  interface.  To  resolve  the  deeper  structures,  a  deep  earth 
model  algorithm  based  upon  standard  seismic  (high  source  angle)  signal 
processing  techniques  is  presented  in  Section  IV  of  this  thesis.  Since  these 
methods  are  industry  standard  in  the  petroleum  exploration  industry,  no 
attempt  will  be  made  to  evaluate  these  techniques.  However,  the  deep  earth 
algorithm  must  be  considered  an  essential  part  of  the  overall  algorithm  for 
Navy  practical  applications. 
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III.  ACOUSTIC  DATA /MODEL  SUMMARY 

A.  EXPERIMENTAL  OVERVIEW  OF  PROJECT  GEMINI 

A  series  of  five  shallow  water,  synthetic  aperture  array  experiments 
collectively  called  Project  GEMINI,  were  collected  off  of  Corpus  Christi, 
Texas  in  the  fall  of  1985.  For  the  purpose  of  this  thesis,  only  those 
measurements  made  at  the  shallow  and  deep  sites  will  be  discussed.  The 
third  site  (Rubano  site)  has  been  studied  extensively  by  previous  efforts  and  a 
well  defined  geoacoustic  characterization  for  the  location  is  described  in 
Matthews  et  al.,  1985  and  Rubano,  1980.  A  recent  analysis  of  the  GEMINI 
region  by  Duarte  (1994)  shows  relatively  good  agreement  between  measured 
and  modeled  transmission  loss  quantities  at  the  Rubano  site.  Given  the 
relatively  benign  environment  of  the  GEMINI  region  and  the  close  proximity 
of  the  array  sites,  it  was  speculated  by  Duarte  (1994)  that  the  geoacoustic 
model,  initially  developed  for  the  Rubano  site,  could  be  successfully  used  at 
the  other  two  sites.  However,  he  concluded  that  the  TL  data  is  particularly 
sensitive  to  the  input  geoacoustic  parameters  and  that  site -specific 
geoacoustic  models  are  required  in  even  in  a  benign,  nearly  range 
independent  environment. 

B.  ENVIRONMENTAL  DESCRIPTION  OF  THE  GEMINI  AREA 

1.  Sound  Speed  Profile 

With  the  exception  of  the  shallow  site,  CTD  and  XBT  casts  were 
made  at  the  commencement  and  completion  of  each  tow  leg.  Two  weeks 
prior  to  the  experiment.  Hurricane  Elena  swept  through  the  Gulf  of  Mexico 
and,  as  a  result,  established  a  well-mixed  water  column  in  the  GEMINI 
region.  Figures  15  and  16  show  the  sound  speed  profiles  for  the  shallow  and 
deep  water  sites,  respectively.  The  sound  speed  profile  for  the  shallow  site 
(Figure  15)  indicates  isothermal  conditions  over  the  entire  water  column. 
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Figure  15.  Sound  speed  profile  for  the  shallow  site  conducted  on  11  September 
1985. 


Figure  16.  Sound  speed  profile  for  the  deep  site  experiment  conducted  on  12 
September  1985. 

The  sound  speed  profile  for  the  deep  site  (Figure  16)  also  shows  the 
presence  of  a  near  isothermal  water  column  down  to  about  41  m  overlying  a 
low  sound  speed  layer  of  approximately  10  m  in  thickness.  This  resxilts  in  a 
thermochne  gradient  of  approximately  -0.7084  8”^  This  negative  sound 
speed  gradient  causes  sound  rays  to  refi’act  downward  and  thus  accentuate 
sound  propagation  in  the  sub-bottom. 
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2.  Bathymetry 

The  Texas  continental  shelf  is  characterized  by  relatively  smooth 
topography  which  gently  dips  toward  the  southeast.  During  Project  GEMINI 
high  resolution  bathymetric  data  were  collected  along  each  of  the  propagation 
paths  as  well  as  at  each  of  the  station  locations.  Bathymetry  versus  range 
from  the  array  locations  are  shown  in  Table  5.  Each  of  the  tow  legs  of  the 
experiment  was  approximately  5000  m  in  length  thus  giving  an  average 
bottom  slope  at  the  shallow  site  of  approximately  0.011°  and  about  0.054°  for 
the  deep  site. 


Bathymetry  at  the  Shallow 
Site 

Bathymetry  at  the  Deep 
Site 

Range  (m) 

Depth  (m) 

Range 

(m) 

Depth  (m) 

0 

21 

0 

62 

543 

21 

767 

61 

2288 

21 

1441 

61 

3783 

20 

2255 

60 

5460 

19 

3425 

60 

4123 

58 

4974 

58 

Table  5.  Bathymetry  for  the  GEMINI  track  experiments  versus  range. 


C.  GEOACOUSTIC  MODELS  FOR  THE  GEMINI  REGION 

1.  Background 

In  order  to  compare  the  measured  TL  to  the  SNAP  and  FEPE  model 
estimates  of  TL,  geoacoustic  models  were  constructed  at  the  shallow  and  deep 
water  sites.  A  geoacoustic  model  describes  the  depth  dependence  of 
compressional  and  shear  wave  speeds,  compressional  and  shear  wave 
attenuation,  and  density  for  each  layer  in  the  seafloor.  In  general,  a 
geoacoustic  model  defines  the  true  thicknesses  and  properties  of  the  sediment 
and  rock  layers  in  the  sea  floor.  Since  a  geoacoustic  model  only  represents 


35 


the  physical  properties  at  a  single  point,  multiple  data  sets  must  be  formed  to 
account  for  the  variations  in  sediment  properties  in  the  horizontal. 

In  a  recent  analysis  of  TL  in  the  GEMINI  area,  Duarte  (1994) 
showed  that  even  in  a  near  range-independent  environment  such  as  the 
GEMINI  area,  model  TL  fluctuations  are  site  specific  and  highly  dependent 
on  the  input  geoacoustic  parameters.  Moreover,  he  found  that  TL  data  and 
model  comparisons  were  in  good  agreement  at  the  Rubano  site  where  a  site 
specific  geoacoustic  model  was  developed  and  incorporated  into  the 
propagation  models.  At  the  shallow  and  deep  sites  the  Rubano  geoacoustic 
data  were  used  in  the  model  calculations  and  rather  poor  agreement  resulted. 
Thus,  it  appears  that  geoacoustic  data  obtained  for  one  location  can  not  be 
used  at  sites  just  12  nm  away,  even  in  a  seemingly  range-independent 
environment.  Because  of  this  horizontal  variability,  it  was  then  necessary  to 
construct  geoacoustic  models  for  the  shallow  and  deep  sites  of  the  GEMINI 
region.  These  are  discussed  in  the  next  section. 

2.  Regional  Geology 

Glaciation  during  the  Pleistocene  era,  0.012-2.8  million  years  before 
present  (mybp),  created  recognizable  imprints  on  Gulf  of  Mexico 
sedimentation  patterns.  Table  6  is  a  generalized  geologic  time  scale  for 
reference.  With  the  growth  and  subsequent  retreat  of  the  ice  sheets,  sea 
level  fell  (regressed)  and  rose  (transgressed)  accordingly.  During  the 
regressive  sequences,  numerous  stream  channels  appeared  along  the 
present-day  Texas  continental  shelf.  Lowering  sea  levels  caused  rapid 
seaward  progradation  of  river  deltas,  which  rapidly  expanded  the  continental 
margins.  During  transgressive  sequences,  when  mean  sea-level  rose,  vast 
areas  of  the  former  coastHne  were  submerged.  The  stream  channels  which 
once  occupied  the  continental  shelf  regions  were  subsequently  filled  with 
fine-grain  sedimentary  deposits. 
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■effSfSHfmSS 

Sea  Level 

■ESSHH 

0.012  mybp  | 

Late  Wisconsin 

Regression 

Mid  Wisconsin 

Transgression 

Early  Wisconsin 

Regression 

Sangamon 

Transgression 

Quaternary 

Late  lUinoian 

Regression 

Pleistocene 

Mid  lUinoian 

Transgression 

Yarmouth 

Transgression 

Kansan 

Regression 

2.1  mybp 

Aftonian 

Transgression 

iiiiiiiiiiiiiiiiiim 

Nebraskan 

Regression 

Pliocene 

BHIH 

■■■■■ 

Table  6.  Geologic  time  scale  for  the  Quaternary  Period  (From  Matthews  et  al., 
1985). 


The  Holocene  age  sediment  layer  represents  the  last  transgressive 
sequence  in  the  post-glacial  period.  This  sequence  is  characterized  by  fine 
grain  silty-clays  and  represents  the  uppermost,  surficial  sediments  in  the 
GEMINI  region.  A  thorough  discussion  of  the  sediment  properties  for  the 
Holocene  sequence  is  presented  in  Matthews  et  al.  (1985).  The  thickness  of 
the  Holocene  layer  varies  substantially  in  the  GEMINI  region.  An  isopach 
chart  of  the  Holocene  age  sediments  is  shown  in  Figure  17  (Berryhill  and 
Tippet,  1981).  In  general,  the  Holocene  layer  thickens  seaward  as  evident  by 
the  seaward  increase  in  two-way  travel  time  from  10  to  40  ms.  The  dashed 
hnes  on  the  chart  represent  the  more  prominent  relict  stream  channels  along 
the  shelf.  The  presence  of  stream  channels  on  multi-channel  seismic  profiles 
is  one  means  of  discerning  the  interface  between  Holocene  age  sediments  and 
the  underlying  late  Wisconsin  sedimentary  sequence  (Matthews  et  al.,  1985). 
Since  these  stream  channels  were  later  filled  by  post-Wisconsin  age 
sediments,  the  Holocene  age  sequences  are  thicker  in  these  localized  regions 
(in  excess  of  12  m)  than  in  the  siirrounding  Holocene  silty-clay  layers. 
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Figure  17.  Sediment  isopach  chart  of  the  GEMINI  region  expressed  in 
milliseconds  of  two  way  travel  time  (broad  solid  line).  Narrow  solid  line 
represents  sediment  deposited  since  the  last  low  stand  of  sea  level.  Dashed  lines 
represent  the  outline  of  localized  thicker  accumulations  of  sediment  in  ancient 
stream  channels.  One  millisecond  is  approximately  equal  to  0.73  m  (Berryhill, 
1981).  GEMINI  array  locations  are  denoted  by  the  ®  symbol. 


The  U.  S.  Naval  Oceanographic  Office  (NAVOCEANO)  in  conjunction 
with  ARL/PSU  collected  16  gravity  cores  in  the  GEMINI  region  (Ross  et  al., 
1978).  The  results  of  their  analysis  indicates  that  the  Holocene  age  silty-clay 
layer  has  a  mean  density  of  1.557  g/cm®  and  an  average  sovmd  speed  ratio 
(i.e.,  ratio  of  sediment  sound  speed  to  bottom  water  sound  speed)  of  0.987. 
These  values  differ  shghtly  from  the  generic  values  presented  in  Hamilton 
(1980)  for  continental  shelf  and  slope  environments.  However,  the  Ross  et  al. 
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data  represents  site -specific  conditions  and  as  such  are  preferred  for  use  in 
the  geoacoustic  models.  Unfortunately,  none  of  the  cores  were  able  to  extend 
deep  enough  to  sample  the  underlying  late -Wisconsin  layers.  Matthews  et  al, 
(1985)  interprets  this  layer  to  be  a  very  fine  sand  or  sandy  silt.  Based  upon 
the  shallow  water  sediment  tables  of  Hamilton  (1980),  the  late-Wisconsin 
sediments  are  assumed  to  have  a  density  of  1.77  g/cm‘^  and  a  sound  speed 
ratio  of  1.080. 

Figure  18  is  a  representative  seismic  profile  (east-west  transect)  for 
the  GEMINI  region.  As  evident  firom  the  figure,  the  seismic  sequence  is 
characterized  by  gently  sloping,  nearly  parallel  reflectors  broken  up  by 
intermittent  growth  faults  and  stream  channels.  The  strong  reflection 
pattern  occurring  at  the  Holocene-late  Wisconsin  interface  is  indicative  of 
the  strong  impedance  contrast  across  the  interface.  Since  acoustic  impedance 
is  the  product  of  a  material’s  density  and  sound  speed,  a  strong  reflection 
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Figure  18.  Multi-channel  seismic  reflection  record  across  the  Texas  continental 
shelf  near  Corpus  Christi.  A  relict  stream  channel  from  the  Late  Wisconsin  trans¬ 
gression  is  shown. 


pattern  on  a  seismic  profile  normally  suggests  a  large  change  in  lithology 
across  an  interface.  The  Holocene  silty-clay  layer  is  characterized  by  a 
relatively  low  sound  speed  and  density  whereas  the  late  Wisconsin  sandy  silt 
layer  is  characterized  by  a  much  greater  sound  speed  and  density  values. 
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This  impedance  contrast  supports  the  well-defined  layering  structure  in 
Figure  18. 

3.  Geoacoustic  Model  Development 

Both  the  shallow  and  deep  water  sites  have  environmental  and 
acoustic  propagation  characteristics  similar  to  the  Rubano  site  which  was 
previously  reported  in  Matthews  et  al.  (1985)  The  sites  differ  only  in  water 
depth  and  in  the  thickness  of  the  Holocene  age  sedimentary  sequence.  Minor 
differences  in  water  column  sound  speed  are  also  apparent,  however  these 
variations  do  not  significantly  impact  the  geoacoustic  characterizations  The 
shallow  site  is  situated  in  about  20.5  m  of  water  and  the  underlying 
Holocene-age  sHty  clay  layer  is  estimated  to  be  8  m  thick.  The  deep  site  is 
situated  in  about  63  m  of  water.  At  the  deep  site  the  Holocene  layer  is 
estimated  to  be  about  24  m  thick. 

At  each  site  the  sound  speed  ratio  for  the  silty-clay  layer  is  assumed 
to  be  0.987,  which  implies  that  the  sound  speed  of  the  silty-clay  layer  is  less 
that  the  overlying  water  column.  Using  methods  described  in  Hamilton 
(1980),  the  compressional  wave  speed  expressed  as  a  function  of  depth  takes 
the  form  of  the  following  regression  formula: 

Vp  =  Vbw  (0.987)  +  1.3  s-RZ)  (46) 

where  Vp  is  the  compressional  wave  speed  in  m/s,  Vbw  is  the  bottom  water 
sound  speed  in  m/s,  and  Z  is  the  depth  of  sedimentary  layer  below  the 
seafloor.  The  product  Vbw(0.987)  represents  the  initial  sound  speed  (Vo)  of  the 
sediment.  The  sediment  so\md  speed  gradient  (1.3  S'^)  was  established  from 
the  NAVOCEANO  cores  and  is  consistent  with  the  near-surface  sound  speed 
gradients  reported  in  Hamilton  (1979).  Because  the  thickness  of  the 
Holocene  sediments  does  not  exceed  a  few  tens  of  meters,  the  higher  order 
terms  of  Equation  (46)  (not  shown)  make  neghgible  contributions  to  the 
sound  speed  profile  in  the  sediment  layer  and  therefore  are  neglected 
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(Matthews  et  al.,  1985).  As  a  result,  the  compressional  wave  speed  in  the 
Holocene  layer  can  be  expressed  as  a  hnear  function  of  depth. 

Thermal  variations  at  the  sea  floor  can  have  a  pronounced 
effect  on  the  upper  sediment  layer  sound  speed  gradient.  If  a  periodic 
thermal  fluctuation  is  applied  at  the  water-sediment  interface,  a  thermal 
wave  will  propagate  into  the  sediment.  The  depth  to  which  the  temperature 
effects  will  be  felt  is  dependent  on  the  sediment  thermal  conductivity  and  the 
heat  flow  (Matthews  et  al.,  1985).  In  a  study  of  the  seasonal  variation  of 
compressional  sound  speed  due  to  temperature  variability  in  the  water 
column,  Raj  an  et  al.  (1992)  showed  that  this  depth  of  influence  was  about  3 
m  below  the  water-sediment  interface.  Therefore,  to  account  for  the  thermal 
variability  at  the  seafloor,  Matthews  et  al.  (1985)  imposed  a  negative  sound 
speed  gradient  in  the  upper  2.5  m  of  the  Holocene  layer.  Below  this  depth 
the  sound  speed  profile  follows  that  expressed  in  Equation  (46). 

The  late-Wisconsin  sediments  can  be  modeled  in  a  similar 
fashion.  Based  upon  the  shallow- water  sediment  tables  of  Hamilton  (1980), 
these  sediments  are  interpreted  to  have  a  sound  speed  ratio  and  density  of 
1.080  and  1.77  g/cm^,  respectively.  In  general,  fine  grained  sands  are 
characterized  by  large  sound  speed  gradients  near  the  surface  decreasing 
sharply  with  depth  (Hamilton,  1980).  As  a  result,  the  sound  speed  profile  as 
a  function  of  depth  is  curve-linear. 

Vp=1681  Z  0.013  (47) 

The  sound  speed  constant  in  Equation  (47)  is  based  on  sediment 
values  at  the  water-sediment  interface.  For  a  buried  layer,  such  as  the  late 
Wisconsin  sand,  compensation  must  be  made  for  the  lithostatic  load  from  the 
overlying  Holocene  sediments.  Initial  values  of  compressional  wave  speed  for 
the  late  Wisconsin  layer  are  estimated  by  matching  the  lithostatic  load 
corresponding  to  the  overlying  Holocene  sediments  to  a  similar  lithostatic 
load  at  depth  in  the  late  Wisconsin  sediment  layer.  At  the  shallow  site  the 
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compressional  wave  speed  at  the  top  of  the  late  Wisconsin  sediment  layer  is 
estimated  to  be  1717.9  m/s.  For  the  deep  site,  with  the  increased  overburden 
pressure,  the  compressional  wave  speed  increases  to  a  value  of  1731.0  m/s. 
Figure  19  illustrates  the  sound  speed  profiles  through  the  different 
lithofacies  for  the  shallow  and  deep  sites.  Of  particular  interest  in  these 
profiles  is  the  relative  thicknesses  of  each  propagating  media.  The  thickness 
of  the  Holocene  sequence  more  than  doubles  in  thickness  between  the 
shallow  and  deep  sites.  As  will  be  shown  later,  this  causes  significant 
fluctuations  in  the  TL  versus  range  curves. 
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Figure  19.  Sound  speed  profiles  at  the  (a)  shallow  and  (b)  deep  sites. 


For  terrigenous  silty-clay  sediments  Hamilton  (1978)  showed 
that  the  sediment  bulk  density  varies  linearly  as  a  function  of  compressional 
wave  speed  for  the  first  few  hundred  meters  below  the  sea  floor.  For  the 
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Holocene  layer,  sediment  density  as  a  function  of  compressional  wave  speed 
was  computed  using  the  equation: 

p=1.135(Vp(Z))- 0.155  (48) 

Here  the  regression  constant  (0.155)  was  adjusted  in  the  equation  to  match 
the  mean  surface  density  of  the  16  cores  collected  by  NAVOCEANO.  For  the 
late  Wisconsin  sequence,  sediment  density  as  a  function  of  compressional 
wave  speed  is  expressed  as: 

p=1.1175(Vp)- 0.1245  (49) 

Since  there  were  no  cores  which  penetrated  into  this  layer,  the 
density  structure  is  speculative  and  is  based  on  historical  trends  from 
Hamilton  (1979). 

Current  estirnates  of  low-frequency  compressional  wave 
attenuation  have  generated  much  debate  in  the  acoustic  community. 
Hamilton  (1972)  reported  the  results  of  in  situ  measurements  of 
compressional  wave  speed  and  attenuation  in  various  sediments  off  San 
Diego.  These  measurements  and  others  from  the  literature  allowed  for 
analyses  of  the  relationships  between  attenuation  and  frequency  and  other 
physical  parameters  to  be  made.  Three  important  relationships  resulted 
from  the  analyses:  (1)  compressional  wave  attenuation  in  marine  sediments 
is  approximately  related  to  the  first  power  of  frequency  in  sands,  muds,  and 
sedimentary  rocks,  (2)  attenuation  can  be  estimated  from  sediment  grain  size 
(porosity)  data,  and  (3)  attenuation  in  marine  sediments  typically  increases 
with  depth  due  to  the  reduction  of  sediment  porosity  to  a  point  where 
overburden  pressure  becomes  the  dominant  effect.  Beyond  this  null  point, 
attenuation  decreases  smoothly  with  depth  and  overburden  pressure 
(Hamilton,  1976). 
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Using  the  methods  proposed  in  Hamilton  (1976),  the  Unear 
relationship  between  compressional  wave  attenuation  and  frequency  takes 
the  form 

a  =  Kpf  "  (50) 

where  a  is  the  attenuation  of  compressional  waves  in  dB/m,  Kp  is  the 
attenuation  in  dB/m-kHz,  f  is  the  frequency  in  kHz,  and  n  is  the  exponent  of 
frequency  (assumed  to  be  1).  The  proportionality  constant  Kp  at  the  water- 
sediment  interface  is  estimated  from  the  sediment  porosity  and  mean  grain 
size.  According  to  Hamilton  (1980),  shallow  water  silty-clay  sediments  have 
an  average  porosity  and  mean  grain  size  (in  phi  units)  of  75.9%  and  8.52(|), 
respectively.  These  values  result  in  an  estimate  of  attenuation  at  the  water- 
sediment  interface  which  is  an  order  of  magnitude  greater  than  pubUshed 
attenuation  results  of  Mitchell  and  Focke  (1980).  Matthews  et  al.  (1985) 
chose  to  use  values  for  sediment  porosity  and  mean  grain  size  between  those 
predicted  for  shallow-water  sediments  by  Hamilton  and  those  values 
indicated  by  Mitchell  and  Focke.  As  a  resvilt,  the  lower  error  bar  estimate  for 
an  8.0(|)  particle  was  used  to  determine  the  proportionaUty  constant.  This 
yielded  an  attenuation  value  of  0.03  dB/m-kHz  at  the  water-sediment 
interface.  An  expression  relating  the  compressional  wave  attenuation  for  the 
Holocene  silty-clay  layer  as  a  function  of  depth  can  thus  be  written  as: 

Kp  =0.030 +  0.0016(Z)  (51) 

The  attenuation  gradient  (0.0016  dB/m^-kHz)  was  estabUshed  in  Matthews  et 
al.  (1985)  and  is  consistent  with  attenuation  gradients  reported  in  Hamilton 
(1976)  for  silty-clays 

Compressional  wave  attenuation  for  the  late  Wisconsin  sand 
layer,  can  also  be  empirically  related  to  grain  size  and  porosity.  In  general, 
sound  attenuation  for  sands  are  much  larger  than  those  of  silts,  clays,  and 
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muds  (Hamilton,  1976).  For  sands,  sound  attenuation  decreases  with  depth 
and  increasing  overburden  pressure.  For  the  late  Wisconsin  sequence,  the 
attenuation  profile  is  calculated  from  the  following  equation  from  Matthews 
et  al.  (1985): 


Kp  =0.29930  -  0.00067(Z)  (52) 

Even  though  the  Holocene  and  late  Wisconsin  age  sediments 
possess  enough  dynamic  rigidity  to  support  shear  wave  propagation,  shear 
wave  speed  and  attenuation  were  not  incorporated  in  the  acoustic  models. 
Given  the  shallow  water  depths,  the  relatively  thin  sediment  cover,  and  the 
low  frequencies  involved  at  the  GEMINI  region,  shear  wave  production  was 
assumed  to  make  a  negligible  contribution  to  the  overall  acoustic  field  (Raj  an, 
1994). 

Geoacoustic  model  tabulations  for  the  shallow  and  deep  water  sites 
are  listed  in  Tables  7  and  8,  respectively.  As  previously  discussed,  the 
models  consist  of  an  upper  sequence  of  Holocene  age  silty-clays  overljdng  a 
late-Wisconsin  age  fine  sand  sequence. 
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Material 

Depth  (m) 

Sound  Speed 

Density 

Sound  Attenuation 

Vp  (m/s) 

(g/cc) 

Kp  (dB/m-kHz) 

Sea  Surface 

0.00 

1531.93 

1.00 

0.000 

1.80 

1531.96 

1.00 

0.000 

2.40 

1539.96 

1.00 

0.000 

3.40 

1545.48 

1.00 

0.000 

4.50 

1545.56 

1.00 

0.000 

5.50 

1545.58 

1.00 

0.000 

6.60 

1545.53 

1.00 

0.000 

7.70 

1545.55 

1.00 

0.000 

8.70 

1545.63 

1.00 

0.000 

Water 

9.40 

1545.65 

1.00 

0.000 

Column 

10.40 

1545.67 

1.00 

0.000 

11.60 

1545.68 

1.00 

0.000 

12.60 

1545.69 

1.00 

0.000 

13.20 

1545.71 

1.00 

0.000 

15.40 

1545.75 

1.00 

0.000 

16.70 

1545.76 

1.00 

0.000 

17.70 

1545.77 

1.00 

0.000 

18.60 

1545.8 

1.00 

0.000 

19.60 

1545.81 

1.00 

0.000 

iVater-Sediment  Interfece 

20.50 

1545.82 

1.00 

0.000 

20.51 

1520.75 

1.57 

0.030 

Holocene  Age 

23.00 

1515.75 

1.56 

0.032 

Silty-Clay 

25.50 

1523.25 

1.57 

0.030 

28.00 

1528.25 

1.58 

0.036 

28.48 

1531.13 

1.58 

0.030 

28.49 

1717.93 

1.79 

0.296 

30.50 

1721.18 

1.80 

0.293 

35.50 

1727.68 

1.81 

0.292 

40.50 

1734.18 

1.81 

0.289 

45.50 

1740.68 

1.82 

0.286 

50.50 

1747.18 

1.83 

0.284 

55.50 

1753.68 

1.84 

0.281 

60.50 

1760.18 

1.84 

0.279 

65.50 

1766.68 

1.85 

0.276 

70.50 

1773.18 

1.86 

0.274 

Late  Wisconsin  Age 

75.50 

1779.68 

1.86 

0.272 

Fine  Sand 

80.50 

1786.18 

1.87 

0.269 

85.50 

1792.69 

1.88 

0.267 

90.50 

1799.18 

1.88 

0.265 

95.50 

1805.68 

1.89 

0.263 

100.50 

1812.18 

1.90 

0.260 

105.50 

1818.68 

1.91 

0.258 

110.50 

1825.18 

1.92 

0.256 

115.50 

1831.68 

1.92 

0.254 

120.50 

1838.18 

1.93 

0.252 

Table  7.  Geoacoustic  Model  for  the  Shallow  Water  Site. 
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Material 

Depth  (m) 

Sound  Speed 
Vp  (m/s) 

Density 

(g/cc) 

Sound  Attenuation 
Kp  (dB/m~kHz) 

Sea  Surface 

0.00 

1544.74 

1.00 

0.000 

1.70 

1544.74 

1.00 

0.000 

2.50 

1544.72 

1.00 

0.000 

3.60 

1544.79 

1.00 

0.000 

4.40 

1544.85 

1.00 

0.000 

5.60 

1544.87 

1.00 

0.000 

6.50 

1544.91 

1.00 

0.000 

7.40 

1544.94 

1.00 

0.000 

8.70 

1544.98 

1.00 

0.000 

9.40 

1545.00 

1.00 

0.000 

10.70 

1545.01 

1.00 

0.000 

11.70 

1545.01 

1.00 

0.000 

12.70 

1544.98 

1.00 

0.000 

13.20 

1544.94 

1.00 

0.000 

15.60 

1544.94 

1.00 

0.000 

16.50 

1544.88 

1.00 

0.000 

17.50 

1544.85 

1.00 

0.000 

18.70 

1544.82 

1.00 

0.000 

19.60 

1544.82 

1.00 

0.000 

20.60 

1544.83 

1.00 

0.000 

21.60 

1544.85 

1.00 

0.000 

22.50 

1544.87 

1.00 

0.000 

23.60 

1544.88 

1.00 

0.000 

25.60 

1544.92 

1.00 

0.000 

26.70 

1544.94 

1.00 

0.000 

27.90 

1544.97 

1.00 

0.000 

Water 

28.40 

1544.97 

1.00 

0.000 

Column 

30.70 

1545.01 

1.00 

0.000 

31.70 

1545.03 

1.00 

0.000 

32.60 

1545.04 

1.00 

0.000 

35.50 

1545.09 

1.00 

0.000 

36.10 

1545.10 

1.00 

0.000 

37.40 

1545.12 

1.00 

0.000 

38.40 

1545.13 

1.00 

0.000 

39.40 

1545.15 

1.00 

0.000 

40.40 

1545.17 

1.00 

0.000 

43.70 

1544.28 

1.00 

0.000 

46.50 

1543.79 

1.00 

0.000 

47.50 

1543.66 

1.00 

0.000 

49.70 

1542.46 

1.00 

0.000 

50.50 

1541.37 

1.00 

0.000 

51.60 

1540.27 

1.00 

0.000 

52.50 

1539.44 

1.00 

0.000 

53.50 

1538.58 

1.00 

0.000 

Table  8.  Geoacoustic  Model  for  the  Deep  Water  Site 
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Material 

Depth  (m) 

Sound  Speed 

Density 

Sound  Attenuation 

Vp  (m/s) 

(g/cc) 

Kp  (dB/m-kHz) 

54.60 

1538.03 

1.00 

0.000 

55.60 

1537.69 

1.00 

0.000 

56.70 

1537.50 

1.00 

0.000 

57.60 

1537.47 

1.00 

0.000 

58.60 

1537.44 

1.00 

0.000 

59.00 

1537.40 

1.00 

0.000 

60.30 

1537.40 

1.00 

0.000 

61.50 

1537.31 

1.00 

0.000 

kVater-Sediment  Interface 

62.34 

1536.99 

1.00 

0.000 

62.35 

1516.03 

1.57 

0.030 

67.34 

1518.53 

1.57 

0.034 

69.84 

1523.53 

1.57 

0.036 

72.34 

1529.03 

1.58 

0.038 

Holocene  Age 

74.34 

1531.63 

1.58 

0.039 

Silty-Clay 

77.34 

1535.53 

1.59 

0.041 

79.84 

1538.78 

1.59 

0.043 

82.34 

1542.03 

1.60 

0.045 

84.84 

1545.28 

1.60 

0.046 

86.00 

1546.97 

1.60 

0.047 

86.01 

1732.07 

1.81 

0.289 

91.00 

1738.60 

1.82 

0.287 

96.00 

1745.13 

1.83 

0.284 

101.00 

1751.66 

1.83 

0.282 

106.00 

1758.20 

1.84 

0.280 

111.00 

1764.73 

1.85 

0.277 

Late  Wisconsin  Age 

116.00 

1771.26 

1.85 

0.275 

Fine  Sand 

121.00 

1777.79 

1.86 

0.273 

126.00 

1784.33 

1.87 

0.270 

131.00 

1790.86 

1.88 

0.268 

141.00 

1803.93 

1.89 

0.263 

151.00 

1816.99 

1.91 

0.259 

156.00 

1823.52 

1.91 

0.256 

161.00 

1830.06 

1.92 

0.254 

Table  8.  Geoacoustic  Model  for  the  Deep  Water  Site  (continued). 


D.  TRANSMISSION  LOSS  MODELS 

In  order  to  qualitatively  assess  the  accuracy  of  the  geoacoustic  models 
(derived  in  the  previous  section),  estimates  of  TL  using  the  geoacoustic 
models  as  input  were  compared  with  the  measured  GEMINI  TL  data.  TL 
must  be  considered  the  measure  of  effectiveness  (MOE)  for  the  performance 
of  the  inverse  technique  developed  in  this  thesis.  If  measured  and  modeled 
TL  agree  exactly  for  two  different  geoacoustic  models,  then  the  differences  in 


the  geoacoustic  parameters  are  considered  acoustically  indistinguishable  for 
a  given  frequency. 

The  normal  mode  model  (SNAP)  and  the  finite  element  parabolic 
equation  model  (FEPE)  were  used  in  this  portion  of  the  analysis  to 
quantitatively  determine  the  impact  of  the  perturbative  technique  to 
estimate  of  the  geoacoustic  model.  The  results  of  these  comparisons  are 
given  in  the  next  section. 

For  the  inverse  theory  problem  developed  in  this  thesis  the  selection  of 
the  proper  acoustic  propagation  model  to  use  is  critical  to  the  success  of  the 
inversion.  The  central  issue  for  this  study  is  the  inversion  of  geoacoustic 
propertied  from  measured  TL  data  assuming  little  to  no  a  priori  geoacoustic 
information  is  available.  Therefore,  it  is  crucial  to  have  an  acoustic  model 
that  is  sensitive  to  even  small  variations  in  geoacoustic  properties  so  that 
when  the  model  and  data  agree  well,  the  geoacoustic  parameters  are 
estimated  accurately. 

Although,  the  GEMINI  region  was  selected  by  the  U.  S.  Navy  by 
virtue  of  its  seemingly  range-independent  environment,  the  geoacoustic 
properties  were  shown  to  be  just  the  opposite.  For  this  reason  SNAP  could 
not  be  used  in  the  inverse  problem  because  of  its  adiabatic  assumption  (i.e., 
not  fully  coupled).  KRAKEN  is  similar  to  SNAP,  except  that  it  is  a  fully 
coupled  normal  mode  model.  KRAKEN  was  used  to  calculate  and  display  the 
normal  mode  functions  and  eigenvalues  (Z(m)  and  km)  that  are  inputs  to  the 
kernel  of  the  Fredholm  integral  in  Equation  (40),  The  FEPE  model  has 
shown  excellent  agreement  in  comparison  with  fully  coupled  normal  mode 
models  and  because  FEPE  is  more  efficient  to  run,  it  was  the  principle  model 
used  to  compare  inverse  results  with  the  measured  TL  data. 
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E.  TRANSMISSION  LOSS  MODEL/DATA  COMPARISONS 

The  GEMINI  experiments  consisted  of  towing  a  narrow-band  CW 
source  at  a  fixed  depth  away  from  a  pair  of  moored  receivers  (Lynch  et  al., 
1991).  An  NRL  J15-3  acoustic  source  with  output  tonals  of  50  and  140.05  Hz 
was  used  during  the  course  of  the  experiment.  The  source  and  receiver 
configuration  is  shown  in  Figure  20  and  in  Table  1.  The  combination  of  the 
two  depth  settings  along  with  the  two  source  frequencies  allowed  for  4 
complex  pressure  field  magnitudes  versus  range  measurements  to  be 
collected  at  each  site.  The  complex  pressure  field  measurements  were 
transformed  into  transmission  loss  (TL)  by  the  relation 

TL(r,z)  =  -201ogl^^  (53) 

|Po| 

where  p(r,z)  is  the  pressure  field  at  range  r  and  po  is  the  acoustic  pressure 
field  at  1  meter  distance  from  the  source.  In  this  section  the  measured  TL 
results  fi:om  the  shallow  and  deep  sites  are  presented  along  with  TL 
estimates  from  FEPE  and  SNAP. 


Figure  20.  Experimental  source  and  receiver  configuration  for  the  GEMINI 
experiments. 
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1.  Shallow  Site 

The  experimental  TL  measurements  for  the  shallow  site  are 
illustrated  in  Figures  21-23  along  with  estimates  of  TL  based  upon  SNAP 
and  FEPE  model  estimates.  The  50  Hz  data  reveal  a  predominantly  two¬ 
mode  interference  pattern  whereas  the  140  Hz  data  show  a  more 
complicated,  multi-mode  interference  pattern.  The  measured  50  Hz  TL  curve 
exhibits  a  high  degree  of  attenuation  attributed  to  a  low  compressional 
velocity  layer  in  the  sub-bottom.  The  curve  also  exhibits  a  high  density  of 
noise  spikes  beyond  about  2500  m  in  range  due  to  local  seismic  profiling 
activity  in  the  region.  These  rapid  fluctuations  should  be  considered  as 
artifacts  in  the  data  and  are  not  expected  to  be  modeled  accurately. 

Figures  21-23  illustrate  the  relatively  close  agreement  between 
the  measured  and  modeled  TL  estimates  using  as  input  the  shallow  site 
geoacoustic  model.  For  the  50  Hz  case  (Figixre  21a),  beyond  2600  m  there  is  a 
substantial  increase  in  TL.  As  previously  shown,  the  sediment  thickness  in 
the  GEMINI  region  is  highly  variable.  The  sharp  increase  in  TL  is  attributed 
to  an  abrupt  increase  in  sediment  thickness  at  about  2600  m  in  range.  This 
feature  is  interpreted  as  one  of  the  numerous  late-Wisconsin  age  stream 
channels  which  criss-cross  the  Texas  continental  shelf  (Matthews  et  al., 
1985).  This  type  of  feature  would  account  for  this  significant,  localized 
increase  in  sediment  thickness.  To  incorporate  this  range-dependency  in  the 
TL  model  estimates,  a  second  geoacoustic  model  was  developed  for  the 
shallow  site.  Based  upon  the  sediment  isopach  chart  (Fig\xre  17),  the 
Holocene  age  silty-clay  layer  was  increased  firom  about  8  m  to  35  m  in 
thickness.  As  evident  fi-om  Figure  23,  this  range-dependent  geoacoustic 
model  improved  the  50  Hz  FEPE  TL  estimates  for  ranges  beyond  2600  m. 
The  140  Hz  signal  (Figure  21b)  was  not  significantly  impacted  by  this  range- 
dependency  since  most  of  the  140  Hz  propagation  occurs  in  the  upper  few 
meters  of  the  sediment. 
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Overall,  the  FEPE  TL  estimates  are  more  accurate  than  the 
SNAP  TL  measurements,  a  result  that  is  in  stark  contrast  to  Duarte’s  (1994) 
conclusions  regarding  TL  model  for  which  he  used  the  Rubano  site 
geoacoustic  data  at  the  shallow  site,  Thus,  it  can  be  concluded  that  the 
shallow  water  TL  model  estimates  are  extremely  sensitive  to  the  geoacoustic 
model  used.  This  conclusion  is  further  illustrated  by  comparing  the  140  Hz 
TL  model  results  from  the  shallow  site  (Figure  21b)  with  those  found  in 
Duarte  (1994)  (Figure  24).  The  lone  difference  between  the  two  FEPE  model 
estimations  rests  with  the  input  geoacoustic  model  used.  As  discussed 
earlier,  the  shallow  site  geoacoustic  model  consists  of  about  8  m  of  near 
surface,  low-velocity  sediments  in  contrast  to  17  m  at  the  Rubano  site 
(Matthews  et  al.,  1985).  The  thicker  low-velocity  sequence  (Rubano  model) 
results  in  more  loss  (>10  dB)  than  that  experienced  by  using  the  shallow  site 
geoacoustic  model.  The  SNAP  TL  estimate  also  shows  reasonably  close 
agreement  with  the  measured  TL  data. 

Figures  21  through  23  suggest  that  the  FEPE  estimates  of  TL  are  in 
closer  agreement  with  the  measured  TL  field  than  those  portrayed  by  SNAP. 
One  possible  reason  for  this  is  that  SNAP  is  not  a  fully  coupled  normal  mode 
model  and,  because  the  geoacoustic  environment  is  range-dependent,  the 
SNAP  model  does  not  permit  energy  to  be  exchanged  between  modes  which 
the  measured  TL  indicates  is  an  important  component. 
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MID-COLUMN  RECEIVER  (50  Hz) 

Range (m) 


SNAP 


FEPE 


GEMINI 


(a) 


MID-COLUMN  RECEIVER  (140  Hz) 


Range  (m) 


Sj^^p  pgpg - GEMINI 


(b) 


Figure  21.  Comparison  of  experimental  and  theoretical  transmission  loss  curves 
for  the  shallow  site  for  a  mid-column  receiver.  Frequency;  50  Hz  (a),  140  Hz  (b). 
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(b) 


Figure  22.  Comparison  of  experimental  and  theoretical  transmission  loss  for  the 
shallow  site  for  a  bottom  receiver.  Frequency:  50  Hz  (a),  140  Hz  (b). 
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Figure  23.  Plot  of  experimental  GEMINI  and  range-dependent  FEPE  model 
transmission  loss  for  the  shallow  site  at  50  Hz. 


Figure  24.  Comparison  of  GEMINI  transmission  loss  data  and  FEPE  transmission 
loss  model  estimates  for  the  shallow  site  based  upon  the  Rubano  site  geoacoustic 
model  as  input  (Duarte,  1994). 
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2.  Deep  Site 

The  deep  site  was  occupied  on  two  different  occasions  during  Project 
GEMINI  (10  and  12  September  1985).  Because  the  data  collected  from  10 
September  are  incomplete  only  the  data  collected  on  12  September  are  used 
in  this  study.  The  experimental  TL  measiorements  for  the  site  are  illustrated 
in  Figures  25-26  along  with  model  estimates  of  TL.  As  with  the  shallow  site 
data,  the  measured  TL  is  characterized  by  a  predominant  two-mode 
interference  pattern  for  the  50  Hz  signal  and  a  multi-mode  interference 
pattern  for  the  140  Hz  signal.  The  FEPE  and  SNAP  model  TL  estimates 
used  the  deep  site  geoacoustic  model  developed  in  the  preceding  section  as 
input.  Unlike  the  shallow  site,  the  propagation  path  at  the  deep  site  was 
treated  as  range-independent. 

Although  there  is  reasonably  good  agreement  in  ampHtude  between 
the  modeled  and  measured  TL  levels  for  the  50  Hz  signal,  both  FEPE  and 
SNAP  fail  to  accurately  model  the  TL  fluctuations  at  the  deep  site.  The  TL 
fluctuations  (Figure  25a)  in  excess  of  20  dB.  These  fluctuations  may  be  due 
to  interfering  normal  modes  at  the  site  (Duarte,  1994). 

Somewhat  better  agreement  is  noted  between  the  modeled  and 
measured  TL  levels  for  the  140  Hz  signal.  The  140  Hz  signal  shows  a  more 
complicated  m\dti-mode  interference  pattern.  As  illustrated  in  Figure  25b, 
TL  fluctuations  are  on  the  order  of  20  dB  in  TL. 
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BOTTOM  RECEIVER  fSO  Hz| 

Range (m| 


0  1000  2000  3000  4000  5000 


0  1000  2000  3000  4000  5000 


SNAP  —  fEPE  -  GEMINI 

(b) 

Figure  26.  Comparison  of  experimental  and  theoretical  transmission  loss  curves 


IV.  PERTURBATIVE  INVERSION  RESULTS 

In  this  section  the  perturbative  inversion  technique  is  applied  to  the 
experimental  data  obtained  during  Project  GEMINI.  As  a  result  of  the 
inversion,  a  new  set  of  geoacoustic  parameters  for  the  shallow  and  deep 
water  sites  are  estimated.  As  discussed  previously,  TL  is  used  as  the 
measure  of  effectiveness  of  the  inverse  geoacoustic  models.  Since  FEPE  was 
shown  to  be  more  accurate  than  SNAP  in  estimating  TL,  only  FEPE  will  be 
used  in  the  comparison  between  measured  and  modeled  TL  values.  The 
major  strength  of  the  inversion  method  is  that  a  priori  information  of  the 
sub-bottom  is  not  required  to  estimate  geoacoustic  parameters.  In  order  to 
demonstrate  this,  we  assume  that  the  geologic  and  geophysical  properties  of 
the  environment  are  only  vaguely  known  in  the  initial  background 
geoacoustic  models. 

A.  SHALLOW  SITE 

The  results  of  Hankel  transforming  the  experimental  pressure  field 
over  the  entire  5000  m  range  at  50  and  140  Hz  are  shown  in  Figure  27.  Since 
the  poles  of  the  Green’s  function  represent  the  eigenvalues  for  discrete  modes 
for  the  particular  waveguide  at  the  given  frequency,  the  location  of  these 
poles  for  the  bottom  receiver  must  be  coincident  with  pole  locations  for  the 
mid-column  receiver.  In  this  section  only  the  plots  for  the  bottom  receiver  are 
illustrated.  However,  in  order  to  form  an  accurate  estimate  of  the  modal 
eigenvalues  of  the  waveguide  (kmeas),  the  results  from  the  two  receivers  are 
averaged  to  form  a  single  estimate  for  each  frequency.  The  eigenvalues  for 
the  two  frequencies  are  tabulated  in  Table  9.  The  50  Hz  plot  (Figure  27a) 
shows  one  distinct  peak  at  a  horizontal  wave  number  of  about  0.19  m“k  A 
less  evident,  very  weak  pole  was  found  to  be  at  a  horizontal  wavenumber  of 
about  0.17  m“k  For  the  140  Hz  case  (Figure  27b),  three  trapped  eigenvalues 
are  readily  identifiable  from  the  plot. 
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Frequency  (Hz) 

Mode  Number 

kmeas  (m  ^) 

50 

1 

0.1896861800 

2 

0.1747074050 

140 

1 

0.5541557518 

2 

0.5431438446 

3 

0.5201824925 

Table  9.  Measured  mode  eigenvalues  for  the  shallow  site. 
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Using  the  modes  shown  in  Table  9,  the  pertubative  inverse  was 
performed  separately  for  each  frequency.  In  order  to  analyze  the  convergence 
of  the  inversion,  two  diverse  initial  background  models  were  used  in  the 
inversion.  The  first  background  model,  BGl,  consisted  of  a  constant  gradient 
{Vc  =  2.0  s"')  background  profile  with  an  initial  sound  speed  of  c(0)  =  1500 
m/s,  and  with  a  constant  density  of  p  =  1.56  g/cm’b  The  second  background 
model,  BG2,  consisted  of  a  constant  gradient  (Vc  =  1.25  s“9  background 
profile  with  an  initial  sound  speed  of  c(0)=1900  m/s  and  with  a  constant 
density  of  p  =  1.56  g/cm'b  Both  background  models  used  the  experimentally 
measured  sound  speed  profile  of  the  water  column  in  the  inversion.  The 
results  of  the  50  Hz  inversion  are  displayed  in  Figure  28.  Also  shown  as 
reference  is  the  geoacoustic  model  developed  in  Section  III. 


Figure  28.  Results  of  the  50  Hz  perturbative  inversion  using  background  models 
BGl  and  BG2.  The  shallow  site  geoacoustic  model  derived  in  Section  III  is  also 
shown. 


Both  profiles  portray  a  lower  sound  speed  sediment  overlying  a  much 
higher  speed  geologic  sequence.  Even  though  there  are  similarities  between 
the  inverse  results  using  BGl  and  BG2,  it  is  apparent  that  the  inverse 
results  of  using  BG2  as  the  background  model  are  in  closer  agreement  with 
the  geoacoustic  model  than  BGl.  This  effect  suggests  that  the  inversion 
process  is  highly  dependent  on  the  initial  constraints  of  the  background 
geoacoustic  model  used  in  the  inversion.  When  this  inversion  algorithm  can 
be  automated,  numerous  initial  estimates  will  be  made  and  the  conversion 
properties  can  be  analyzed  in  more  detail.  If  only  a  few  initial  estimates  are 
made,  the  low  eigenvalues  can  be  used  to  assess  stabihty  of  the  inversion  of 
the  spectral  decomposition  in  the  Backus-Gilbert  solution.  The  higher  the 
magnitude  of  the  low  level  eigenvalues,  the  more  stable  and  accurate  the 
solution. 

Using  the  same  initial  background  geoacoustic  models,  the  results  of 
recovering  the  sound  speed  profile  for  the  140  Hz  signal  are  shown  in  Figure 
29.  As  the  number  of  modes  and  frequency  increases,  large  fluctuation  in  c(z) 
result.  To  remedy  this  situation,  a  simple  “box-car”  average  filter  was 
apphed  to  c(z).  As  with  the  50  Hz  case,  the  140  Hz  inversion  suggests  a  two 
layer  geologic  sequence  consisting  of  a  lower  speed  layer  overlying  a  higher 
speed  layer. 
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Figure  29.  Results  of  the  140  Hz  perturbative  inversion  using  background  models 
BGl  and  BG2.  The  shallow  site  geoacoustic  model  derived  in  Section  III  is  also 
shown. 

It  is  important  to  emphasize  that  the  inversion  process  did  not 
incorporate  any  a  priori  information  of  the  geologic  and  geophysical 
parameters  of  the  seafloor.  Geophysical  parameters  such  as  sediment 
density  and  compressional  attenuation  were  held  constant  throughout  the 
process.  In  practice  however,  certain  knowledge  of  the  sediment  properties 
may  be  obtained  and  thus  make  the  problem  less  undetermined.  For 
example,  by  incorporating  knowledge  of  the  sound  speed  ratio,  the  profiles  in 
Figures  28  and  29  can  be  adjusted  such  that  the  sound  speed  profile 
intersects  the  seafloor  at  the  measured  sound  speed  of  the  sediment. 
Through  conventional  seismic  methods,  information  on  the  deep  layering 
structure  may  be  established  which  will  place  lower  bounds  on  the 
perturbation  integral.  Although  not  used  in  these  analyses  since  only  low 
grazing  angle  data  were  collected,  the  deep  earth  algorithm  described  in  a 
subsequent  section  could  have  accurately  identified  the  sedimentary  layers 
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shown  in  Figure  18.  By  incorporating  this  type  of  information  into  the 
inversion,  the  solution  is  less  undetermined  and  converges  rapidly  to  a 
solution. 

In  order  to  test  the  accuracy  of  the  shallow  site  inversion,  model 
estimates  of  TL  were  compared  with  measured  TL  values.  For  the  50  Hz 
case,  the  results  from  BGl  were  used  as  the  input  to  FEPE.  It  is  evident 
from  Figure  28  that  the  derived  sound  speed  profile  suggests  a  two  layer 
sedimentary  sequence  consisting  of  a  lower  speed  sediment  overlying  a 
higher  speed  layer.  The  rather  sharp  transition  in  sound  speed  at  about  35 
m  in  depth  indicates  a  rather  abrupt  change  in  lithology  which  we  interpret 
as  the  change  from  the  sdty-clay  layer  to  the  sand.  Therefore,  generic  values 
of  sediment  density  and  compressional  attenuation  were  incorporated  in  the 
inverse  geoacoustic  model  to  account  for  this  transition.  For  the  140  Hz  case 
(Figure  29),  the  results  from  BGl  were  modified  in  a  similar  fashion  to 
incorporate  the  change  in  density  and  attenuation  with  depth. 

Figures  30  and  31  show  a  comparison  of  the  TL  curves  measured  at 
the  shallow  site  with  the  TL  curves  found  through  pertubative  inversion 
technique  for  the  same  receiver  combinations  as  in  the  previous  section.  Also 
shown  in  the  figures  are  TL  versus  range  results  using  the  geoacoustic  model 
derived  in  Section  II.  As  evident  from  the  curves,  the  perturbative  inversion 
model  results  in  estimates  of  TL  that  are  as  good  or  in  some  case  better  than 
using  a  high  resolution  geoacoustic  model  as  input  to  the  acoustic 
propagation  models. 
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Figure  30.  Plot  of  GEMINI  TL  data  versus  FEPE  using  as  input  the  perturbative 
inverse  solution  and  a  site  speciJQc  geoacoustic  model  for  the  mid-column  receiver, 
shallow  site.  Frequency  :  50  Hz  (a),  140  Hz  (b). 
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Figure  31.  Plot  of  GEMINI  TL  data  versus  FEPE  using  as  input  the  perturbative 
inverse  solution  and  the  site  specific  geoacoustic  model  for  the  bottom  receiver, 
shallow  site.  Frequency  :  50  Hz  (a),  140  Hz  (b). 
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B.  DEEP  SITE 


The  Green’s  function  versus  horizontal  wavenumber  plots  for  the  deep 
site  are  shown  in  Figures  32  and  33.  As  with  the  shallow  site,  these  plots  are 
the  result  of  Hankel  transforming  the  experimental  pressure  field  over  the 
entire  propagation  path.  Typically,  as  frequency  and  water  depth  increase, 
the  number  of  trapped  mode  substantially  increases.  Tbis  is  made  apparent 
by  the  large  number  of  poles  of  the  Green’s  function.  The  eigenvalues  for  the 
two  frequencies  are  tabulated  in  Table  10.  The  50  Hz  plot  (Figure  32)  shows 
two  distinct  peaks  at  horizontal  wavenumbers  of  about  0.203  m‘i  and  0.197 
m'k  A  weak  pole  is  also  evident  at  about  a  horizontal  wavenumber  of  about 
0.190  m~’.  The  140  Hz  plot  shows  four  distinct  modes  situated  at  horizontal 
wavenumbers  of  about  0.57,  0.55,  0.54,  and  0.52  m“k  However,  only  the 
latter  three,  more  energetic  modes  were  actually  used  in  the  inversion,  since 
the  pole  at  the  horizontal  wavenumber  of  0.57  did  not  correlate  with  results 
from  the  mid-column  receiver. 
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Figure  32.  Magnitude  of  the  depth  dependent  Green’s  function  at  50  Hz  for  the 
deep  site  experiment.  Only  the  bottom  receiver  is  shown. 
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Green's  Function  versus  Horizontal 
Wavenumber 
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Figure  33.  Magnitude  of  the  depth  dependent  Green’s  function  at  140  Hz  for  the 
deep  site  experiment.  Only  the  bottom  receiver  is  shown. 


Frequency  (Hz) 

Mode  Number 

kmeas  (m 

50 

1 

0.20322170900693 

2 

0.19715011547344 

3 

0.19012702078522 

140 

1 

0.55415575176589 

2 

0.54314384460141 

3 

0.52018249243189 

Table  10.  Measured  mode  eigenvalues  for  the  deep  site. 


Using  the  modes  shown  in  Table  10,  the  perturbative  inverse  was 
performed  separately  for  each  frequency  at  the  deep  site.  Unlike  the  shallow 
site,  the  inversion  process  at  the  deep  site  was  extremely  sensitive  to  the 
initial  background  model  used.  The  initial  background  models  BGl  and  BG2 
used  in  the  previous  section  resulted  in  an  unstable  inverse  solution.  One 
explanation  for  this  is  that  low  order  eigenvalues  in  the  Backus-Gilbert 
decomposition  result  in  an  undetermined  inverse  matrix  which  causes  the 
solution  to  become  unstable.  However,  several  experimental  iterative  runs 
showed  that  a  background  model  starting  with  a  reduced  initial  sound  speed 
of  c(0)  =  1400  m/s,  and  a  constant  gradient  of  Vc  =  2.0  s“^  coupled  with  a 
constant  density  of  p  =  1.56  g/cm'^  provided  a  stable  sound  speed  profile.  The 
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results  of  the  50  Hz  inversion  are  displayed  in  Figure  34  along  with  the  site 
specific  geoacoustic  model  developed  previously. 


Figure  34.  Results  of  the  50  Hz  perturbative  inversion  at  the  deep  site.  The  deep 
site  geoacoustic  model  derived  in  Section  III  is  also  shown. 


Using  the  same  initial  background  geoacoustic  model  as  for  50  Hz,  the 
results  of  the  140  Hz  perturbative  inversion  are  shown  in  Figure  35.  In  order 
to  stabilize  the  fluctuations  in  c(z),  a  “box-car”  average  was  once  again 
applied  to  the  140  Hz  inverse  solution.  Both  the  50  Hz  and  140  Hz  sound 
speed  profides  suggest  that  the  upper  sedimentary  layer  at  the  deep  site 
extends  deeper  than  indicated  in  the  geoacoustic  model.  This  is  a  valid 
result,  since  as  discussed  previously,  the  thickness  of  the  Holocene  age  silty- 
clay  layer  varies  substantially  in  the  GEMINI  region. 
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Figure  35.  Results  of  the  140  Hz  perturbative  inversion  at  the  deep  site.  The  deep 
site  geoacoustic  model  derived  in  Section  III  is  also  shown. 

Using  values  consistent  with  Hamilton  (1978,  1980),  the  above  sound 
speed  profiles  were  adjusted  to  account  for  varying  sediment  density  and 
compressional  attenuation  with  depth.  The  resulting  geoacoustic  models 
w'ere  used  as  input  to  FEPE.  Figures  36  and  37  show  a  comparison  of  the  TL 
curves  measured  at  the  deep  site  with  FEPE  estimates  of  TL  using  the 
results  of  the  corrected  inversion  as  input.  As  evident  from  the  figures,  the 
estimates  of  TL  generated  by  the  perturbative  inverse  solution  are  extremely 
consistent  with  the  measured  TL  values. 
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Figure  3®.  Comparison  of  GEMINI  TL  data  and  FEPE  TL  mode!  estimates  using  as 
input  the  perturbative  inverse  solution  and  the  site  specific  geoacoustic  model  for 
the  deep  site  for  a  mid-column  receiver.  Frequency:  50  Hz  (a),  140  Hz  (b). 
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Figure  37«  Comparison  of  GEMINI  TL  data  and  FEPE  TL  model  estimates  using  as 
input  the  perturbative  inverse  solution  and  the  site  specific  geoacoustic  model  for 
the  deep  site  for  a  bottom  receiver.  Frequency:  50  Hz  (a),  140  Hz  (b). 
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C.  DEEP  EARTH  SEISMIC  MODEL 

The  inverse  technique  developed  in  this  thesis  relies  principally  on 
shallow  grazing  angle  (<  20°  )  propagation  in  estimating  the  near  surface 
geoacoustic  parameters.  Typically,  using  this  inversion  method,  geoacoustic 
parameters  can  be  accurately  resolved  to  about  4/3  of  an  acoustic  wavelength 
down  from  the  seafloor.  Inasmuch  as  low  angle  propagation  paths  are 
important  for  this  inversion  technique  and  more  importantly  in  typical  ASW 
scenarios,  the  technique  is  not  effective  in  identifying  underlying  sediment 
layers  beyond  4/3  acoustic  wavelengths. 

To  augment  the  capabihty  of  the  inverse  technique,  conventional 
seismic  reflection  methods  are  proposed  to  resolve  these  deeper  sediment 
structures.  These  methods  typically  rely  on  higher  propagation  grazing 
angles  (  >  45  °)  and  thus  are  able  to  penetrate  deeper  into  the  seafloor. 
Therefore,  a  deep  earth  model  is  required  to  account  for  these  deeper  geologic 
horizons.  A  deep  earth  model  was  not  applied  to  experimental  data  in  this 
thesis;  however  a  brief  discussion  is  provided  in  this  section  to  show  how  this 
information  may  be  incorporated  into  the  perturbative  inverse  technique. 

The  proposed  deep  earth  model  is  developed  from  high  angle  (>  45  °) 
propagation  paths  and  low  (<  25  Hz)  source  frequencies.  The  supporting  data 
can  be  collected  simultaneously  as  the  data  collected  for  the  inverse 
technique  by  collecting  an  adequate  sample  of  short  range  (i.e.,  high  angle) 
time-series.  Subsequently,  the  processing  then  follows  standard  seismic 
signal  processing  procedures  as  follows: 

Step  1:  The  time  series  data  are  stacked  or  gathered  using  common 
midpoint  (CMP)  stacking  algorithms. 

Step  2:  The  data  set  is  deconvolved  to  remove  artifacts  from  the  data 
such  as  water  borne  reverberation. 

Step  3:  Normal  move  out  (NMO)  corrections  are  appHed  to  the  CMP 
stacks. 
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Step  4:  Frequency  filter  the  data  with  a  time  dependent,  bandpass 
filter.  The  bandpass  filters  are  constructed  in  a  manner  to  allow  lower 
frequencies  to  pass  as  two  way  travel  time  increases. 

Step  5.  In  range -dependent  sub-bottom  environments,  time  and  depth 
migration  schemes  should  be  employed  to  identify  lateral  discontinui¬ 
ties  in  the  strata  (i.e.,  faults,  intrusions,  etc.) 

Step  6;  Display  data  in  a  distance  versus  two-way  travel  time  plot 
such  that  geologic  Hthologies  and  structures  can  be  identified. 

If  steps  1  through  6  were  applied  to  the  GEMINI  data,  the  sediment 
interface  between  the  Holocene  age  silty-clay  layer  and  the  late  Wisconsin 
age  sand  layer  could  be  easily  discerned  by  source  firequencies  as  high  as  50 
Hz.  Incorporating  the  exact  position  of  the  layer  into  the  inversion  would 
reduce  the  undetermined  state  of  the  inverse  solution  and  ultimately  would 
improve  the  final  product.  The  Navy  now  has  the  low  firequency  capability 
(e.g.,  HLF-1  sonar  on  submarines)  on  board  and  thus  the  deep  earth  model  is 
very  practical  and  cost  effective.  The  deep  earth  model  should  be  considered 
part  of  the  inverse  technique  algorithm  recommended  for  use  in  the  fleet. 


76 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  technical  approach  in  this  thesis  addressed  two  methods  of  estimat¬ 
ing  geoacoustic  parameters  for  acoustic  propagation  modeling.  In  the  first 
part  of  the  study  propagation  model  estimates  of  TL  were  compared  with 
measured  low-frequency  acoustic  propagation  data  collected  during  Project 
GEMINI.  Location-specific,  high  resolution  geoacoustic  models  of  the  sub¬ 
bottom  were  input  to  the  acoustic  models  in  the  first  part  of  the  study.  The 
geoacoustic  models  derived  in  this  section  were  constructed  in  manner 
described  in  Hamilton  (1980)  and  made  use  of  available  geologic  and 
geophysical  data  The  acoustic  models  examined  in  this  portion  of  the  study 
were  the  SACLANTCEN  normal  mode  acoustic  propagation  model  (SNAP), 
and  the  NRL  finite  element  parabohc  equation  model  (FEPE).  The 
comparisons  were  made  over  two  source  frequencies  (50  and  140  Hz)  at  two 
depth  settings  (mid-column  and  near-bottom  setting).  Overall,  it  was  shown 
that  FEPE,  using  as  input  the  detailed  site  specific  geoacoustic 
characterizations,  resulted  in  estimates  of  TL  that  were  more  accurate  than 
those  calculated  by  SNAP  and  were  reasonably  consistent  with  the  measured 
TL  field.  Moreover,  it  was  shown  that  TL  estimates  are  extremely  sensitive 
to  the  input  environmental  parameters.  Even  in  a  seemingly  range - 
independent  acoustic  environment  such  as  off  the  Texas  coast,  a  small 
change  in  a  geoacoustic  parameter  (e.g.,  sediment  thickness)  can  result  in  a 
significant  change  in  TL. 

In  the  second  part  of  this  study,  the  perturbative  inversion  technique 
(Rajan,  1987)  was  used  to  estimate  the  geoacoustic  parameters  for  the  two 
GEMINI  sites.  The  objective  of  perturbative  inversion  is  to  estimate  the 
depth  variation  in  geoacoustic  parameters  in  the  sub-bottom  (i.e.,  sound 
speed,  density,  attenuation,  etc.)  from  the  differences  between  the  horizontal 
wave  numbers  from  measured  data  and  the  corresponding  wavenumbers 
computed  for  a  background  geoacoustic  model.  To  illustrate  the  efficacy  of 
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the  inverse  method,  the  initial  background  geoacoustic  model  used  in  the 
inversion  assumed  no  a  priori  information  of  the  seafloor.  The  resultant 
inverse  solution  was  used  to  predict  the  geoacoustic  properties  at  each  of  the 
sites.  As  a  measime  of  effectiveness,  estimates  of  TL  using  the  inverse 
results  (modeled  with  FEPE)  were  directly  compared  to  the  measured 
GEMINI  TL  values.  The  final  results  were  in  excellent  agreement  with  the 
measured  TL  and  the  resulting  TL  estimates  were  as  good  as  or  better  than 
the  TL  estimates  obtained  from  detailed,  site-specific  geoacoustic  models. 

By  assuming  no  a  priori  information  in  the  inversion  process  and 
subsequently  obtaining  such  excellent  results  in  TL  addresses  an  extremely 
critical  problem  for  the  U.  S.  Navy  in  shallow  water  regions.  In  many  of  the 
strategic,  littoral  environments  knowledge  of  the  geoacoustic  properties  of 
the  sub-bottom  is  extremely  limited.  Direct  geologic  and  geophysical  surveys 
in  these  regions  to  determine  relevant  geoacoustic  environmental  parameters 
are  both  costly  and  time  consuming  and  generally  not  consistent  with  the 
rapid  response  requirements  of  today’s  fleet.  Therefore,  the  use  of  perturba¬ 
tive  inversion  technique  is  an  accurate,  cost  effective  and  practical 
alternative  to  obtaining  acoustic  performance  prediction  estimates  in  these 
regions.  Data  necessary  for  the  inversion  can  be  easily  collected  using 
existing  Navy  platforms  and  hardware. 

Automation  of  the  inversion  algorithm  presented  in  this  thesis  could 
allow  for  an  entire  range  of  initial  background  models  to  be  implemented  in 
the  inversion  process.  At  present,  one  background  model  is  assumed  and 
after  each  iteration  of  the  inversion,  a  set  of  eigenvalues  and  eigenfunctions 
are  generated  by  the  normal  mode  program  (KRAKEN).  However,  in  a  fully 
automated  mode,  an  entire  range  of  background  models  with  differing  initial 
sound  speeds  and  sound  speed  gradients  could  be  analyzed  by  the  inversion 
routine  consecutively.  By  applying  this  method,  a  family  of  sound  speed 
profiles  can  be  formulated.  This  is  especially  applicable  in  high  noise 
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environments  where  the  eigenvalues  of  the  inverse  matrix  are  exceptionally 
small. 

Currently,  there  are  numerous  existing  pressure  field  data  sets  (e.g., 
Bearing  Stake,  NATIVE  I,  Southwest  Florida  Sea  Test,  etc.)  that  the 
inversion  method,  as  described  in  this  thesis,  can  be  readily  applied.  To  test 
the  robustness  of  the  inversion,  several  environmental  scenarios  should  be 
examined  ranging  from  a  weakly  range-dependent  environment  such  as  the 
Southwest  Florida  shelf  to  a  more  geologically  complex  region  such  as  the 
Northern  Arabian  Sea  (Bearing  Stake  exercise).  Finally,  it  is  recommended 
that  an  at-sea  experiment  be  conducted  using  fleet  assets  to  collect  the 
required  pressure  field  data. 
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